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ABSTRACT 


Vincent  (2006)  developed  a  technique  for  remotely  measuring  Aerosol  Optical  Depth 
(AOD)  using  commercial  high-resolution  satellite  imagery.  This  technique  measured  the 
radiance  difference  between  a  building  shadow  and  an  adjacent  sunlit  region  with  the 
same  surface  reflectance  to  calculate  Total  Optical  Depth  (TOD).  AOD  is  then 
determined  by  subtracting  Rayleigh  scattering  from  TOD.  The  procedure  for  making  this 
calculation  was  time  consuming,  particularly  locating  suitable  shadows  within  the  region 
of  interest.  This  paper  outlines  a  fully  automated  method  of  performing  the  AOD 
calculation  and  examining  shadow  properties.  The  automated  method  relies  on  high- 
resolution  Digital  Surface  Model  (DSM)  data  collected  using  a  Light  Detection  and 
Ranging  (LIDAR)  sensor  coupled  with  sun  and  satellite  geometry  to  identify  shadow 
regions.  Configuration  settings  allowed  for  specific  regions  in  the  shadow  and  sunlit  area 
to  be  selected  before  determining  their  respective  radiances.  Finally,  a  technique  for 
aligning  the  satellite  and  DSM  pixels  was  developed  to  correct  for  small  differences 
between  the  datasets.  Results  from  the  automated  method  were  compared  with 
AERONET  data  for  validation.  The  automated  method  using  WorldView-1  and 
QuickBird  imagery  worked  best  at  Solar  Village,  Saudi  Arabia  and  an  area  northeast  of 
Washington,  D.C.,  which  included  the  Goddard  Space  Flight  Center.  Testing  of 
IKONOS  multispectral  imagery  suggested  the  resolution  is  inadequate  in  urban  settings. 
Testing  in  areas  that  included  downtown  regions  in  Houston,  TX  and  Baltimore,  MD 
identified  weaknesses  in  the  alignment  algorithm. 
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I.  INTRODUCTION 


Aerosols  play  an  important  part  in  the  earth’s  radiation  budget  directly  (through 
scattering  and  absorption  of  incoming  solar  radiation  and  outgoing  longwave  radiation) 
and  indirectly  (by  affecting  cloud  properties)  (IPCC  2007).  The  overall  impact  depends 
on  the  aerosol  optical  properties  that  vary  with  wavelength,  relative  humidity,  and  the 
horizontal  and  vertical  distribution,  all  of  these  varying  with  time  (Hayward  and  Boucher 
2000,  IPCC  2001).  Although  aerosols  can  contribute  positively  to  global  heating,  their 
overall  effect  is  a  net  cooling  effect.  This  cooling  counteracts  25-50%  of  the  heating 
from  anthropogenic  greenhouse  gases  (IPCC  2007).  Similarly,  aerosols  significantly 
attenuate  transmission  of  electromagnetic  radiation  in  the  visible  and  near-infrared  (NIR) 
through  scattering  and  absorption,  limiting  the  effectiveness  of  intelligence  collection  and 
weapon  systems  that  operate  in  those  portions  of  the  spectrum  by  reducing  the  amount  of 
radiation  reaching  the  sensors.  Aerosol  optical  depth  (AOD)  is  a  measure  of  the  opacity 
of  the  atmosphere  due  to  aerosols.  It  is  wavelength-dependent  and  is  directly 
proportional  to  the  size  and  concentration  of  the  aerosol  particles  present.  Accurately 
characterizing  the  radiative  forcing  due  to  aerosols  is  important  when  running  climate 
models  (IPCC  2007).  However,  the  significant  spatial  and  temporal  variability  makes  it 
difficult  to  directly  measure  aerosol  properties  (size  and  concentration)  in  the  atmosphere 
globally.  As  a  result,  numerous  algorithms  using  satellite  imagery  to  determine  aerosol 
properties  have  been  developed. 

Most  satellite  techniques  for  characterizing  aerosols  with  visible  and  NIR 
wavelength  sensors  rely  on  surfaces  with  low  reflectance  in  the  visible  channels  or  over 
surfaces  where  detailed  knowledge  of  the  surface  reflectance  exists.  A  few  satellites 
(Polarization  and  Directionality  of  the  Earth’s  Reflectance  -  POLDER  and  Multi-angle 
Imaging  Spectro-Radiometer  -  MISR)  have  dedicated  sensors  for  measuring  aerosol 
properties  (IPCC  2007).  These  sensors  use  multiple  look  angles  to  model  the  atmosphere 
but  are  limited  to  low  surface  reflectance  regions  (Veefkind  et  al.  1998,  Martonchik  et  al. 
2004,  Deuze  et  al.  2000).  Vincent  (2006)  developed  a  technique  to  characterize  aerosols 
using  commercial  high-resolution  satellite  data  in  the  visible  and  NIR  portions  of  the 
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spectrum,  over  any  terrain,  without  prior  knowledge  of  or  constraints  on  surface 
reflectance.  This  technique,  henceforth  referred  to  as  the  shadow  technique  or  method, 
compares  the  radiance  inside  the  shadow  with  radiance  just  outside  the  shadow  (sunlit)  to 
calculate  AOD.  Although  promising,  Vincent’s  method  was  labor  intensive,  making 
large  case  studies  and  practical  application  difficult.  Since  that  time,  additional  research 
has  refined  the  technique  and  sped  up  the  process  (Rivenbark  2009,  Sweat  2008, 
Dombrock  2007,  Evans  2007).  However,  those  efforts  stopped  well  short  of  full 
automation. 

This  research  develops  and  tests  a  method  called  the  Fully  Automated  Shadow 
Technique,  or  FASTEC,  of  fully  automating  the  shadow  technique  for  any  sensor  or 
wavelength  provided  a  high  resolution  surface  elevation  map  is  available.  In  addition  to 
identifying  and  calculating  AOD  for  specific  shadows  or  every  usable  shadow  in  an 
image,  FASTEC  enables  detailed  examination  of  shadow  and  sunlit  region  characteristics 
for  determining  pixel  selection  criteria.  Chapter  II  describes  the  shadow  technique  as 
developed  by  Vincent  (2006).  It  also  lists  and  briefly  describes  other  methods  of 
measuring  AOD.  Chapter  III  outlines  the  automated  approach  and  its  implementation. 
Chapter  IV  details  the  results  of  testing  the  automated  method.  Chapter  V  discusses 
conclusions  and  offers  suggestions  for  future  research  with  the  automated  method. 
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II.  BACKGROUND 


A.  OTHER  SATELLITE-BASED  AOD  RETRIEVAL  METHODS 


Numerous  satellite-based  techniques  for  measuring  AOD  exist.  These  techniques 
typically  work  with  specific  wavelength  bands,  minimum  horizontal  resolutions,  and 
surface  reflectance  ranges.  Figure  1  is  a  global  seasonal  average  composite  of  AOD  from 
the  Moderate  Resolution  Imaging  Spectrometer  (MODIS)  sensor,  demonstrating  the  high 
spatial  and  seasonal  variability  of  AOD  globally.  Additionally,  Figure  1  shows  large 
areas  where  AOD  is  not  calculated  because  existing  MODIS  techniques  do  not  work  over 
desert  or  snow-covered  regions  due  to  high  surface  reflectance. 


Figure  1. 
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Aerosol  optical  depth  at  0.55  /urn  using  the  MODIS  sensor  (from  IPCC  2007). 
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Surface  reflectance  is  negligible  over  water  for  near-infrared  (NIR)  wavelengths 
and  AOD  is  calculated  using  the  dark  object  method,  which  is  described  in  Chapter  II.A.2 
(von  Hoyningen-Huene  2003).  Here  total  optical  depth  (TOD)  is  computed  by  treating 
the  low  surface  reflectance  as  negligible  and  assuming  all  radiance  received  by  the 
satellite  is  scattered  spaceward  by  aerosol  particles  and  gases.  After  removing  the 
contribution  from  Rayleigh  scattering  (depends  on  surface  pressure),  the  contribution 
from  aerosol  scattering  dominates  the  remaining  signal.  Over  land,  determining  optical 
depth  is  less  straightforward  since  the  contribution  from  surface  reflectance  is  non-zero 
and  varies  significantly  with  surface  type  and  wavelength.  To  calculate  AOD  over  land, 
more  complex  methods  like  contrast  reduction,  dark  object,  and  multi-angle  were 
developed,  and  these  are  described  as  follows: 

1.  Contrast  Reduction  Method 

The  contrast  reduction  method  computes  optical  depth  from  the  measured 
radiance  difference  between  a  distinct  surface  and  its  surrounding  area  (Odell  and 
Weinman  1975).  It  requires  prior  knowledge  of  the  two  surface  albedos,  sun  and 
satellite  zenith  angles,  and  lookup  tables  based  on  the  mean  surface  albedo,  the  viewing 
zenith  angle,  and  a  contrast  transmission  function  to  compute  optical  depth.  As  defined 
by  Odell  and  Weinman  (1975),  the  contrast  transmission  function  is  proportional  to  the 
difference  between  the  reflected  radiance  off  the  surfaces  and  inversely  proportional  to 
the  difference  between  their  surface  albedos.  Odell  and  Weinman  (1975)  showed  that  it 
is  possible  to  exploit  shadows  cast  by  large  objects  to  detennine  optical  depth  without 
prior  knowledge  of  the  surface  features.  This  forms  the  basis  of  the  shadow  method 
developed  by  Vincent  (2006).  Later,  Kaufman  and  Joseph  (1982)  applied  a  similar 
approach  for  regions  where  the  surface  reflectance  undergoes  a  step-wise  change,  such  as 
along  a  shoreline  or  forest  border,  to  compute  optical  depth  and  surface  reflectivity.  This 
eliminates  the  previous  requirement  that  surface  reflectance  be  known.  Additionally, 
using  two  different  surface  regions  allows  this  approach  to  be  independent  of  satellite 
calibration  and  produces  more  accurate  results  than  earlier  contrast  reduction  methods. 
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2. 


Dark  Object  Method 


The  dark  object  method  is  usable  for  wavelengths  with  low  surface  reflectance. 
In  such  locations  it  is  possible  to  estimate  the  small  portion  of  the  top-of-atmosphere 
(TOA)  signal,  which  is  due  to  surface  reflection.  Kaufman  and  Sendra  (1988)  detail  a 
rapid,  automated  algorithm  to  compute  optical  thickness,  using  the  dark  object  method 
over  regions  covered  by  dark  vegetation.  Their  approach  assumes  the  surface  reflectance 
of  those  regions  to  estimate  optical  thickness  and  was  accurate  to  within  0.20  for  band  1 
(0.5-0. 6  pm)  and  0.05  for  band  2  (0.6-0. 7  pm).  The  dark  object  method  is  ineffective 
over  surfaces  with  high  surface  reflectance,  such  as  snow  or  sand.  Other  dark  object 
algorithms,  such  as  that  described  by  von  Hoyningen-Huene  et  al.  (2003),  employ  the 
same  fundamental  idea. 

3.  Multi-angle  method 

Multi-angle  methods  use  special  sensors  that  examine  a  region  from  multiple  look 
angles.  Sensors  such  as  POLDER,  MISR,  or  Along  Track  Scanning  Radiometer  2 
(ATSR-2)  characterize  aerosols  by  using  the  different  look  angles  to  measure  surface 
reflectance  and  use  lookup  tables  to  model  the  atmosphere  (Veefkind  et  al.  1998, 
Martonchik  et  al.  2004,  Deuze  et  al.  2000).  For  example,  the  dual  view  algorithm, 
developed  by  Veefkind  et  al.  (1998),  uses  two  look  angles  to  separate  the  atmospheric 
and  surface  contributions  to  the  satellite  radiance.  The  multi-angle  method,  though  quite 
accurate,  includes  a  step  in  the  algorithm  to  identify  and  ignore  regions  with  high  surface 
reflectance  (Robles  Gonzalez  et  al.  1997). 

B.  THE  SHADOW  METHOD 

Except  where  otherwise  cited,  the  information  provided  in  this  description  of  the 
shadow  method  was  taken  from  Vincent  (2006).  The  shadow  method  exploits  high- 
resolution  imagery  to  measure  the  radiance  difference  between  building  shadows  and 
adjacent  sunlit  regions  to  calculate  total  optical  depth.  Several  commercially-available 
satellites  have  sufficient  horizontal  resolution,  from  46  cm  for  the  WorldView-1 
panchromatic  channel  to  2.44  m  for  the  QuickBird  multispectral  channels  (DigitalGlobe, 
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cited  June  2010a/b),  to  use  the  shadow  method.  Figure  2  depicts  the  basic  components  of 
radiance  received  at  a  satellite,  from  a  shadow  and  an  adjacent  sunlit  region  where  both 
occur  over  an  identical  surface.  Direct  transmission,  diffuse  transmission,  and  diffuse 
reflection  components  make  up  the  sunlit  signal,  while  the  shadow  signal  lacks  the  direct 
transmission  component.  For  a  homogeneous  surface,  where  a  portion  is  covered  by 
shadow,  subtracting  the  shadow  radiance  from  the  sunlit  radiance  isolates  the  direct 
transmission  component  of  the  signal,  the  portion  of  the  signal  that  is  the  key  to 
calculating  optical  depth. 


Figure  2.  Depiction  of  the  difference  between  the  shadow  and  sunlit  signal  received  by 
a  high-resolution  satellite.  Subtracting  the  two  signals  isolates  the  direct 
transmission  component  of  the  signal,  which  can  be  used  to  calculate  total  optical 

depth  (from  Vincent  2006). 

1.  General  Calculation 

The  equation  for  calculating  TOD  is  derived  from  Liou’s  (2002)  development  of 
the  principles-of-invariance  (POI)  method  to  describe  the  multi-scattering  problem. 
Vincent  (2006)  provides  a  detailed  derivation  of  the  governing  equation: 
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where  50  is  total  optical  depth,  p0  is  the  cosine  of  the  solar  zenith  angle,  p  is  the  cosine  of 


the  sensor  zenith  angle,  rs  is  the  surface  reflectance,  his  the  mean  aerosol  reflectance 
(MAR),  F0  is  the  spectral  solar  irradiance  for  the  specific  channel,  and  Ld  is  the  difference 
between  the  sunlit  and  shadow  radiance.  As  the  radiance  difference  increases,  TOD 
decreases.  Every  term  except  Ld  in  Equation  (1)  is  either  included  with  the  image  (in  a 
metadata  file)  or  can  be  estimated.  Y  and  rs  were  estimated  using  the  following 
equations: 
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where  d>  is  the  single  scatter  albedo,  P  is  the  scattering  phase  function  (represented  by 
the  Henyey-Greenstein  phase  function),  g  is  the  asymmetry  parameter  (0.65),  ©is  the 
scattering  angle,  and  (p',^')and  (p",^")are  the  cosine  of  the  incident  and  resultant  zenith 


and  azimuth  angles  (Henyey  and  Greenstein  1941).  The  single  scatter  albedo  value 
varies  regionally,  but  Vincent  used  values  of  0.88  for  Beijing,  China  and  0.95  for  other 
locations.  Additionally,  Vincent  set  the  asymmetry  parameter  value  to  0.65  to  account 
for  both  aerosol  (0.6-0. 7)  and  Rayleigh  particles  (~  0). 


Y  approximates  the  diffuse  sky  radiance.  Although  the  diffuse  field  consists  of 
aerosol  and  molecular  affects,  the  molecular  affects  were  ignored  when  establishing  the 
asymmetry  parameter  in  Equation  (3)  because  aerosol  affects  dominate  at  all  but  the 
lowest  optical  depths.  Had  the  asymmetry  parameter  included  the  molecular  affects,  its 
value  would  decrease.  Vincent  concluded  this  would  cause  the  scattering  phase  function 
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and  MAR  to  increase,  thereby  lowering  the  AOD  and  increasing  the  low  bias  of  the 
shadow  method.  Therefore,  molecular  affects  were  only  accounted  for  during  the  final 
partitioning  of  the  total  optical  depth. 

To  estimate  AOD,  the  contribution  from  Rayleigh  scattering,  or  molecular 
Rayleigh  optical  depth  (5r),  must  be  calculated  and  subtracted  from  the  TOD.  The  effect 
of  molecular  absorption  is  considered  negligible  since  most  satellite  sensor  wavelength 
bands  operate  in  atmospheric  transmission  windows.  The  molecules  responsible  for 
Rayleigh  scattering  are  standard  atmospheric  constituents  and  may  be  calculated  globally 
for  a  given  wavelength  using  equations  defined  by  Russell  et  al.  (1993): 

5r  (A,)  =  (0.00864  +  6.5  x10~6h)  Ub(x)  (p  /  p0)  (5) 

b(k)  =  3.916  +  0.074k  +  0.050/k  (6) 

where  H  is  the  height  above  sea  level  of  the  radiometer  in  kilometers,  p  is  the 
atmospheric  pressure  at  H  in  hPa,  and  po  is  the  sea-level  reference  pressure  (1013.25 
hPa).  The  depolarization  factor  (0.00864)  was  chosen  based  on  the  best  fit  data  (Frolich 
and  Shaw  1980,  Young  1980).  Using  Equations  (5)  and  (6),  5r  was  approximated  using 
the  center  effective  wavelength  of  the  satellite  channel. 

2.  Sensitivities 

Vincent’s  sensitivity  analysis  of  key  terms  in  Equation  (1)  demonstrated  the 
shadow  method  is  most  sensitive  to  the  radiance  difference  between  sunlit  and  shadow 
regions  and  surface  reflectance  and  much  less  sensitive  to  7,  solar  zenith  angle,  and 
sensor  zenith  angle.  Additionally,  Vincent  detailed  the  expected  impact  on  the  shadow 
method  when  cumuliform,  stratiform,  or  cirrus  clouds  were  present.  In  the  case  of 
cumuliform  or  stratifonn  clouds,  images  with  cloud  coverage  of  more  than  20%  were 
avoided.  Thin  cirrus  was  more  difficult  to  detect  using  visible  channels  and  resulted  in 
retrieved  radiances  at  the  sensor  being  higher,  causing  an  overall  high  bias  to  the 
calculated  AOD. 
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3. 


Limitations  and  Target  Selection  Criteria 


Vincent  also  found  several  limitations  to  the  shadow  technique.  First,  the 
minimum  allowed  radiance  difference  was  10  W  m"  sr'  nm'  (RU)  and  the  minimum 
allowed  surface  reflectance  was  0.15.  These  thresholds  were  set  because  of  the  high 
sensitivity  to  changes  in  radiance  difference  and  surface  reflectance  below  those  values. 
Such  sensitivity  leads  to  large  errors  in  AOD  calculations  with  small  errors  in  radiance 
difference  or  surface  reflectance.  The  range  of  calculable  AOD  was  set  from  0.1  to  2.0. 
The  minimum  value  was  established  to  avoid  negative  AOD  after  subtracting  5r 
(Rivenbark  2009).  The  maximum  was  detennined  by  calculating  TOD  using  the 
combination  of  allowable  values  that  gave  the  maximum  result.  Furthermore,  when 
optical  depth  is  above  2.0  the  surface  is  obscured  and  shadows  are  indistinguishable. 

The  method  also  required  careful  target  selection  in  urban  environments  and  only 
blue  and  NIR  channels  provided  reliable  results  over  grass  backgrounds.  This  resulted  in 
several  recommended  target  selection  criteria.  These  criteria  included  using  large 
shadows  over  a  homogeneous  surface  with  a  minimum  shadow  size  of  4  pixels  wide  and 
4  pixels  deep  (from  about  1  to  32  square  meters,  depending  on  satellite  resolution),  using 
shadows  that  are  at  least  1  km  away  from  other  structures,  and  avoiding  highly  reflective 
shadowing  structures  and  the  subsequent  adjacency  effects.  Vincent  suggested  two 
approaches  to  automating  the  technique,  dynamic  and  fixed.  The  dynamic  approach 
relied  on  isolating  shadows  of  opportunity  while  the  fixed  approach  required  the 
preidentification  of  buildings  meeting  specified  criteria  for  analysis. 

4.  Other  Observations 

Regardless  of  channel  or  background,  the  shadow  technique  was  a  manually 

intensive  procedure  with  more  time  being  spent  identifying  good  shadows  than 

performing  the  actual  calculation.  Process  streamlining  perfonned  by  Rivenbark  (2009) 

and  Dombrock  (2007)  took  steps  to  improve  the  calculation  time  while  exploring  aspects 

of  shadows  that  contributed  to  error  in  optical  depth.  In  each  study,  errors  in  the  code 

written  in  Interactive  Data  Language  (IDL)  impacted  their  findings.  For  example,  when 

calculating  the  scattering  phase  function  as  part  of  the  Y  calculation  in  Equation  (3),  a 
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syntax  error  in  code  raised  the  denominator  to  the  power  of  1  instead  of  1.5.  The  overall 
impact  of  this  mistake  is  difficult  to  determine  mathematically,  since  the  error  can  cause 
either  an  increase  or  a  decrease  in  the  value  of  the  scattering  phase  function,  depending 
on  the  scattering  angle.  Also,  additional  iterations  in  the  azimuth  angle  loop  nests  in 
Equation  (2)  contributed  to  the  low  bias  by  causing  higher  values  of  7  and  lower  values 
of  optical  depth. 

C.  SHADOW  DETECTION  IN  SATELLITE  IMAGERY 


The  key  step  in  developing  a  fully  automated  approach  to  the  shadow  method  is 
finding  suitable  shadows  within  the  image.  One  satellite-based  method  of  finding 
shadows  uses  color  space  indices  (Sarabandi  et  al.  2004,  Gevers  and  Smeulders  1999). 
The  color  space  indices  are  defined  as  follows: 

_  f  R 

C,  =  arctan  - 

vmax(G,B) 

_  f  G  ) 

C,  =  arctan  - 

vmax(R,B)y 


C,  =  arctan  - 

^  max(R,  G) 

where  R,  G,  and  B  represent  the  radiance  from  the  red,  green  and  blue  color  bands, 
respectively.  Coupling  the  output  from  these  indices  with  a  3x3  texture  filter  clearly 
outlines  the  boundary  between  a  shadow  and  sunlit  region.  Sarabandi  et  al.  (2004) 
determined  that  C3  worked  best  for  locating  shadows.  Due  to  more  Rayleigh  scattering 
of  blue  light  than  red  or  green  light,  C3  is  larger  in  shadow  regions  than  sunlit  regions. 

Although  capable  of  identifying  shadows,  this  satellite -based  approach  is  poorly 
suited  for  automating  the  shadow  method  for  two  reasons.  First,  C3  may  be  less  effective 
in  regions  with  large  aerosol  concentrations  since  the  amount  of  scattered  red  and  green 
light  will  increase.  Second,  this  technique  lacks  necessary  infonnation  to  select  sunlit 
pixels  with  the  same  surface  reflectance  as  shadow  pixels,  particularly  when  using 
panchromatic  channels,  where  surface  reflectance  estimation  techniques  are  not  available. 
Similar  problems  would  likely  arise  with  other  satellite-based  approaches.  The  next 
section  describes  the  FASTEC  procedures  for  identifying  usable  shadows  and  calculating 
AOD. 
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III.  METHODOLOGY 


An  automated  approach,  based  on  trigonometry,  is  created  using  the  shadow 
technique  and  overcoming  the  shortfalls  of  satellite-based  shadow  identification 
techniques  mentioned  in  Chapter  II,  Section  C.  Knowing  pixel  elevations  and  the 
azimuth  and  elevation  of  the  sun  and  satellite,  the  process  for  identifying  shadow  and 
sunlit  pixels  is  relatively  straightforward.  If  the  line-of-sight  between  the  sun  and  a  pixel 
(called  the  target  pixel)  is  blocked  by  a  higher  elevation  pixel,  a  target  pixel  is  in  shadow. 
If  not,  a  target  pixel  is  sunlit  (see  Figure  4).  Similarly,  the  new  approach  ensures  the 
satellite  does  not  collect  radiance  off  an  obstruction,  instead  of  the  target  by  ensuring  the 
line-of-sight  between  the  satellite  and  the  target  pixel  is  not  blocked. 

A.  TARGET  SELECTION  CRITERIA  TESTING 

Past  research  identified  several  target  selection  criteria  for  identifying  good 
shadow  candidates  for  the  shadow  method  (Vincent  2006,  Dombrock  2007,  Rivenbark 
2009).  Criteria  include  avoiding  partial  shadow  pixels  along  the  shadow  edge  and 
avoiding  pixels  where  surrounding  objects  limit  diffuse  sky  radiance.  Detailed 
descriptions  of  how  the  automated  approach  addresses  these  concerns  are  described  in 
Chapter  III,  Section  C.  The  hypotheses  listed  below  describe  key  target  selection  criteria. 

Before  these  criteria  were  tested,  variables  used  to  discriminate  between  shadow 
generators  were  tested  to  eliminate  tree  shadows  from  the  sample.  Tree  shadows  may 
allow  some  sunlight  to  penetrate  through  the  leaf  canopy  into  the  shadow,  increasing  the 
shadow  radiance,  invalidating  the  theory  behind  the  shadow  technique  (see  Figure  2).  The 
first  of  these  variables  sets  the  minimum  height  of  the  shadow  generator  and  the  second 
sets  the  minimum  number  of  pixels  required  for  AOD  calculation.  To  detennine  the 
combination  of  these  variables  that  limits  the  sample  to  buildings,  different  values  of  the 
minimum  shadow  generator  height  were  applied  to  the  GSFC  tile  and  the  output 
examined  to  find  the  value  that  leaves  only  building  shadows.  Furthermore,  plots  of 
AOD  for  the  entire  DSM  were  created  to  show  how  it  varies  as  tree  shadows  are 
introduced  to  the  sample. 
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Figure  3.  Example  of  pixel  rejections  from  FASTEC.  Orange  represents  the  building, 

black  is  the  shadow,  and  gray  is  rejected  shadow.  Panel  1  is  the  full  shadow,  2  is 
after  edge  rejection,  and  3  is  after  building  proximity  rejection. 

1.  Hypothesis  I:  It  Is  Important  to  Avoid  Partial  Shadow  Pixels  Along 
the  Shadow  Edge 

Subjective  examination  of  satellite  imagery  shows  that  radiance  values  transition 
across  several  pixels  along  the  shadow  edge.  In  the  transition  region,  shadow  radiance  is 
slightly  higher  and  sunlit  radiance  is  slightly  lower  than  pixels  that  are  fully  in  shadow  or 
sunlit.  Including  these  pixels  in  the  optical  depth  calculation,  lowers  the  radiance 
difference  and  increases  AOD.  The  edge  pixel  rejection  variable  enables  the  automated 
technique  to  eliminate  partial  shadow  pixels  along  the  shadow  edge  (see  Figure  3).  To 
test  this  hypothesis,  plots  of  radiance  based  on  distance  from  shadow  edge  were  reviewed 
to  determine  the  extent  of  the  influence  of  the  shadow  edge  on  the  shadow  region. 
Additionally,  the  edge  pixel  rejection  threshold  was  varied  and  the  resulting  AOD  values 
compared  for  individual  shadows  and  an  entire  DSM  file.  Chapter  III,  Section  B.5.C. 
describes  the  algorithm  FASTEC  employs  to  eliminate  partial  shadow  edge  pixels. 

2.  Hypothesis  II:  It  Is  Important  to  Choose  Shadow  Pixels  Far  from  the 
Shadow  Generator  to  Calculate  Optical  Depth 

Based  on  the  past  research,  it  is  suggested  the  shadow  region  should  only  include 
pixels  far  from  the  shadow  generator.  In  theory,  reflection  off  surrounding  objects  may 
increase  pixel  radiance  while  blocking  by  surrounding  objects  may  decrease  pixel 
radiance.  Regardless  of  which  modification  dominates,  the  total  signal  is  no  longer 
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composed  of  the  basic  components  listed  in  Figure  2.  By  rejecting  a  percentage  of  the 
shadow  pixels  near  the  shadow  generator  based  on  the  size  of  the  overall  shadow  (see 
Figure  3)  this  adjacency  affect  is  reduced.  To  test  this  hypothesis,  the  correlation 
between  the  pixel  radiance  and  distance  from  the  shadow  generator  was  calculated  for 
several  shadows.  Also,  the  building  percentage  rejection  threshold  was  varied  on 
individual  shadows  and  on  an  entire  DSM  file  to  examine  the  impact  on  AOD.  Chapter 
III,  Section  B.5.a.  describes  the  algorithm  FASTEC  employs  to  eliminate  pixels  too  close 
to  the  shadow  generator. 

3.  Hypothesis  III:  Sky  Dome  Blockage  Is  Important  When  Choosing 
Valid  Shadow  and  Sunlit  Pixels  for  Calculating  Optical  Depth 

Based  on  the  same  reasons  as  described  for  Hypothesis  II,  selecting  pixels  based 
on  sky  dome  blockage  may  be  important  to  the  shadow  technique.  In  order  to  determine 
the  influence  of  sky  dome  blockage  on  the  radiance  of  a  pixel,  the  correlation  between 
them  was  calculated  for  several  shadows.  Additionally,  AOD  was  calculated  for  different 
maximum  sky  dome  blockage  thresholds  for  individual  shadows  and  an  entire  DSM  file. 
Chapter  III,  Section  B.5.b.  describes  the  algorithm  FASTEC  employs  to  eliminate  pixels 
based  on  sky  dome  blockage. 

4.  Hypothesis  IV:  Eliminating  the  Extreme  Radiance  Values  from  the 
Sample  Improves  AOD  Calculations 

Radiance  variation  exists  in  the  shadow  and  sunlit  regions  due  to  variations  in 
surface  reflectance  and  position  differences  between  the  satellite  and  DSM  files.  Surface 
reflectance  variations  may  be  caused  by  different  surfaces  or  different  objects  such  as 
trees,  sidewalks,  or  cars.  Position  differences  between  the  satellite  and  DSM  files  are 
partially  corrected  by  the  offset  correction  algorithm,  but  some  error  in  position,  rotation, 
or  distortion  may  still  exist.  Therefore,  FASTEC  eliminates  percentages  of  the  extreme 
radiance  values  within  shadow  and  sunlit  regions  to  try  and  ensure  the  sampled  regions 
are  representative  of  the  true  shadow/sunlit  radiances.  Assuming  FASTEC  correctly 
identifies  the  shadow,  the  bulk  of  the  pixels  within  the  samples  are  representative  of  the 
true  shadow/sunlit  radiances.  Sorted  values  of  pixel  radiance  were  plotted  to  demonstrate 
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this  problem.  Also,  thresholds  for  eliminating  extreme  values  of  radiance  from  both  the 
shadow  and  sunlit  region  were  varied  with  resulting  AOD  values  plotted  for  both 
individual  shadows  and  across  an  entire  DSM.  Chapter  III,  Section  B.5.d.  describes  the 
method  used  to  eliminate  extreme  radiance  pixels. 

5.  Hypothesis  V:  It  Is  Possible  to  Get  an  Idea  of  Uncertainty  When 
Calculating  AOD  on  Entire  Scenes. 

After  setting  target  selection  criteria  based  on  the  tests  described  by  Hypotheses  I- 
IV,  calculations  on  an  entire  scene  will  provide  a  number  of  shadow/sunlit  pairs  and 
computed  AOD  values.  The  range  of  AOD  values  in  the  scene  will  give  an  idea  of  the 
uncertainty  of  the  automated  technique.  To  test,  AOD  values  were  calculated  for  all 
available  imagery  using  entire  DSM  files.  In  cases  where  ground  truth  data  is  also 
available,  AOD  values  for  each  shadow  were  plotted  against  ground  truth  and  bias  was 
calculated.  The  data  is  analyzed  by  region  (for  each  DSM  file)  as  well  as  across  all 
regions  (all  DSM  files).  In  the  case  of  multispectral  imagery,  plots  of  AOD  between 
channels  were  examined.  Since  the  overlap  between  most  of  the  available  imagery  and 
the  DSM  files  contains  several  qualified  shadow  generators,  a  range  of  AOD  values  was 
available  for  most  test  cases. 

6.  Hypothesis  VI:  The  DSM  Can  Have  Significantly  Lower  Resolution 
Than  the  Satellite  Data  and  Still  Work  Effectively 

The  computational  expense  of  the  automated  algorithm  is  most  influenced  by  the 
size  of  the  DSM  file,  limiting  the  size  of  the  region  over  which  the  automated  method  can 
be  run  operationally.  Using  lower  resolution  imagery  allows  calculations  over  a  larger 
area  without  increasing  computation  time.  For  example,  increasing  the  horizontal 
resolution  of  the  DSM  file  from  one  meter  to  two  meters,  quadruples  the  areal  coverage 
of  the  DSM  with  the  same  computational  cost.  As  long  as  the  resolution  is  sufficient  to 
identify  shadow  regions  with  several  pixels,  it  is  believed  the  resulting  AOD  values  will 
be  comparable  to  the  full  resolution  DSM.  To  test  this  hypothesis,  results  with  varying 
horizontal  resolution  DSM  files  were  compared  to  the  full  resolution  results. 
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B. 


TRIGONOMETRIC  SHADOW  DETECTION  METHOD 


Figure  4.  Trigonometric  approach  employed  by  FASTEC  to  find  shadow  pixels  based 
on  elevation.  Red  arrows  indicate  shadow  pixels  whose  path  to  the  sun  is  blocked 
by  a  higher  elevation  pixel.  The  green  arrow  indicates  an  unblocked  sunlit  pixel. 
Black  pixels  are  obstructed  from  view  of  the  satellite  by  the  building. 


A  high-resolution  surface  elevation  map,  known  as  a  Digital  Surface  Model 
(DSM),  provides  the  critical  data  for  FASTEC.  Using  this  data  instead  of  satellite  data, 
ensures  that  the  automated  approach  works  in  the  presence  of  obscurations  such  as  clouds 
and  aerosols.  To  execute,  FASTEC  requires  overlapping  DSM  and  satellite  files  as 
inputs.  Output  includes  a  list  of  shadows  with  total  optical  depth  (TOD)  and  AOD  values 
for  each  qualifying  shadow  in  the  DSM  file  and  general  details  of  shadow/sunlit  selection 
criteria.  Figure  5  depicts  a  top-level  flowchart  of  the  automation  program.  Subroutines 
for  each  block  in  this  flowchart  are  described  in  detail  in  the  rest  of  this  section. 
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Figure  5.  Top-level  flowchart  of  FASTEC. 

1.  Preparation 

The  preparation  stage  is  made  up  of  two  steps  with  step  one  being  performed  by 
the  user  prior  to  running  the  program  and  step  two  being  performed  by  the  automation 
program  during  execution.  Before  the  automated  method  can  begin,  the  DSM  and 
satellite  files  are  converted  into  flat  file  format,  composed  of  a  header  file  and  a  data  file. 
This  makes  it  easy  to  retrieve  key  header  information.  The  details  of  this  step  as  well  as  a 
description  of  the  input  files  are  detailed  in  Chapter  III,  Section  C.  1  and  2. 

The  automated  portion  of  the  preparation  stage  begins  by  reading  in  a  host  of 
variables  from  a  configuration  file.  The  configuration  file  allows  the  user  to  quickly 
change  variable  values  and  eliminates  the  need  to  recompile  the  program  whenever  these 
variable  values  are  changed.  For  this  program,  the  configuration  file  variables  include  the 
input  and  output  file  names,  shadow  and  sunlit  modification  variables,  offset  correction 
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and  dome  blockage  criteria,  and  AOD  calculation  variables  (see  Appendix  A  for  a 
description  of  each  variable).  Next,  key  information  from  the  DSM  and  satellite  flat  file 
headers  are  read  into  the  program  including:  the  number  of  rows  and  columns  of  data, 
the  number  of  channels  of  data,  the  data  storage  format,  the  latitude  and  longitude  of  the 
top  left  pixel,  and  pixel  size  in  both  the  north-south  and  west-east  directions. 
Additionally,  necessary  infonnation  from  the  image  metadata  file  is  read  including:  the 
absolute  calibration  factor  and  effective  bandwidth  of  each  channel,  the  date  and  time  of 
the  image,  the  satellite  name,  and  the  sun  and  satellite  azimuths  and  elevations.  Equation 
(8)  converts  Digital  Numbers  (DN)  into  spectral  radiance  using  absolute  calibration 
factor  and  effective  bandwidth: 


L,(i)  =  DN, 


absCalFactj 

effBandwidtlr 


(8) 


where  Lx(i)is  the  spectral  radiance,  i  is  the  channel  number,  and  absCalFactj  and 


effBandwidth;  are  the  absolute  calibration  factor  and  effective  bandwidth,  respectively. 
DN  is  essentially  the  relative  spectral  radiance  and  is  based  on  the  amount  of  energy 
detected  by  the  sensor.  Additionally,  the  values  for  §r  and  spectral  solar  irradiance  are 
chosen  from  Table  1  during  this  step.  6R  values  were  calculated  by  taking  a  response- 
weighted  average  using  Equations  (5),  (6),  and  (9): 


8rt  — 


j*8R(A.)0(A.)dA. 


(9) 


where  5RT  is  the  total  MROD  across  the  channel  and  (p(k)  is  the  specific  spectral 


response  for  that  wavelength  and  satellite  channel.  An  example  of  the  spectral  response 
curve  for  the  QuickBird  satellite  is  provided  in  Figure  6.  The  summation  in  Equation  (9) 
was  performed  using  a  50  mn  wavelength  interval  across  the  channel.  Spectral  solar 
irradiance  was  calculated  using  the  same  approach  where  the  wavelength-specific 
spectral  solar  irradiances  were  taken  from  Wehrli  (1985). 
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Sat 

Band 

Wavelength 
Range  (nm) 

Center  Effective 
Wavelength  (nm) 

Molecular  Rayleigh 
Optical  Depth 

Spectral  Solar 
Irradiance 
(W  m'2  nm'1) 

Quickbird 

Blue 

450-520 

486 

0.172 

1940 

Green 

520-600 

547 

0.105 

1851 

Red 

630-690 

652 

0.051 

1593 

Near-Infrared 

760-900 

804 

0.026 

1127 

Panchromatic 

445-900 

689 

0.061 

1558 

Ikonos 

Blue 

445-51 6 

497 

0.165 

1916 

Green 

506-595 

560 

0.097 

1835 

Red 

632-698 

666 

0.047 

1539 

Near-Infrared 

757-853 

792 

0.027 

1170 

Panchromatic 

526-929 

727 

0.048 

1468 

WV1 

Panchromatic 

400-900 

667 

0.070 

1587 

Table  1.  List  of  Molecular  Rayleigh  Optical  Depth  and  Spectral  irradiance 

constants  by  satellite  and  channel. 


The  last  steps  in  the  preparation  stage  expand  the  row  and  column  sub-region  (if  a 
sub-region  was  defined)  and  reads  in  the  entire  DSM  file.  Expanding  the  sub-region 
enables  FASTEC  to  identify  usable  sunlit  pixels  outside  the  user-defined  area  to  pair  with 
shadows  inside  the  sub-region.  The  entire  DSM  file  data  is  stored  in  an  array  for  the  step 
described  in  the  next  section. 


- Pan 

- Blue 

—  Green 

-  .Red 

-  -  NIR 


Wavelength 


Figure  6.  Relative  spectral  response  for  the  QuickBird  satellite  derived  from  data 

provided  by  DigitalGlobe. 
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2.  Identify  Shadow  and  Sunlit  Pixels 


This  step  identifies  shadow  and  sunlit  pixels  and  determines  which  pixels  are 
blocked  from  the  view  of  the  satellite.  The  results  are  stored  in  an  array  and  used  by 
subsequent  steps  in  the  program.  The  results  array  contains  one  of  three  possible  values: 
distance  from  the  shadow  generator  for  in-shadow  pixels,  zero  for  sunlit  pixels,  or  invalid 
for  pixels  blocked  from  view  by  the  satellite. 


Figure  7.  Calculation  of  in-shadow  pixel  h(rtar,Ctar).  The  height  of  the  beam  (hb) 
between  the  target  pixel  and  the  sun  is  less  than  the  height  of  pixel  h(rios,cios).  Diff 
is  the  difference  between  h(rios,cios)  and  hb. 


Pixels  are  considered  shadow  whenever  the  line-of-sight  between  the  target  pixel 
and  the  sun  is  blocked  by  a  higher  elevation  pixel  (see  Figure  7).  The  elevation  of  the 
line-of-sight  between  the  sun  and  target  pixel  is  computed  using  the  below  equations: 


los  ^tar 


dist  =  ^[dy(rlos  -rtar)]'  +[dx(c 

h  b  =  h  (  bar  > ' Car  )  +  diSt  X  taI1  (  eleVsun  ) 


(10) 

(11) 


where  hb  is  the  elevation  at  pixel  (rios,  cios)  of  the  line-of-sight  between  the  sun  and  the 
target  pixel  (rtar,Ctar),  dx  (dy)  is  the  pixel  size  in  meters  in  the  east/west  (north/south) 
direction,  h(rtar,ctar)  is  the  height  of  the  target  pixel,  and  elevsun  is  the  elevation  angle  of 
the  sun  in  radians.  This  calculation  is  repeated  for  every  pixel  on  the  azimuth  between 
the  target  pixel  and  the  sun  until  hb  exceeds  the  elevation  of  the  tallest  pixel  in  the  DSM. 
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Shadow  pixels  store  the  distance  from  the  target  to  the  shadow  generator  in  the 
result  array.  Some  shadow  pixels,  such  as  those  having  two  tall  buildings  between  them 
and  the  sun,  may  have  multiple  shadow  generators.  In  that  situation,  the  result  array 
value  is  determined  by  the  pixel  with  the  largest  difference  (diff)  between  the  shadow 
generator  height  h(rios,cios)  and  hb  (see  Figure  7).  This  ensures  the  stored  value  coincides 
with  the  shadow  generator  that  places  the  target  pixel  deepest  in  shadow.  This  approach 
will  cause  an  occasional  good  shadow  pixel  to  be  marked  invalid  when  eliminating 
unwanted  pixels  (see  Chapter  III,  Section  B.5.a).  Any  pixel  whose  line-of-sight  is 
unblocked  is  considered  sunlit  and  is  labeled  with  a  zero  in  the  result  array. 

After  identifying  each  pixel  in  the  DSM  as  shadow  or  sunlit,  the  view  between  the 
satellite  and  the  ground  is  considered.  This  ensures  that  pixels  representing  objects 
obstructing  the  surface  when  viewed  from  the  satellite  are  identified  and  labeled  invalid 
in  the  result  array.  The  procedure  for  perfonning  this  test  uses  a  similar  approach  as 
described  above  with  the  satellite  location  substituted  for  that  of  the  sun.  Aside  from 
that,  the  only  difference  is  that  any  obstructed  surface  region  is  immediately  marked 
invalid  in  the  result  array.  When  this  step  completes,  the  result  array  contains  pixels 
identified  by  distance  from  shadow  generator  for  shadow  pixels,  zero  for  sunlit  pixels,  or 
invalid  for  pixels  obstructed  from  view  of  the  satellite. 

3.  Identify  Coherent  Shadows 

Once  shadow  and  sunlit  pixels  are  identified,  FASTEC  groups  adjacent  shadow 
pixels  as  one  shadow  and  assigns  that  shadow  a  number,  sequentially  from  one,  for 
identification.  The  procedure  begins  by  examining  the  result  array  from  the  previous 
section  (any  pixel  whose  value  is  greater  than  zero).  Once  a  shadow  pixel  is  found  a 
recursive  subroutine  is  called  to  find  neighboring  shadow  pixels.  Any  neighboring 
shadow  pixel  that  is  found  is  identified  as  belonging  to  the  same  shadow  and  the  same 
recursive  subroutine  is  repeated  for  that  shadow  pixel.  This  process  is  repeated  until 
every  shadow  pixel  has  been  assigned  a  shadow  number.  To  improve  perfonnance  in 
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later  steps,  the  range  of  rows  and  columns  and  the  pixel  count  for  each  coherent  shadow 
are  also  stored.  The  coherent  shadow  array  contains  positive  integers  for  shadow  pixels, 
zero  for  sunlit  pixels,  and  invalid  for  other  pixels. 

4.  Correct  for  Pixel  Offset  Between  DSM  and  Satellite  File 

The  latitude  and  longitude  of  each  point  in  the  DSM  and  satellite  file  must  match 
precisely  for  FASTEC  to  properly  measure  shadow  and  sunlit  radiance  and  calculate 
AOD.  Due  to  the  poor  geolocation  of  the  satellite  imagery  used  in  this  study,  pixel 
locations  within  these  files  rarely  match  and  correction  is  required.  To  demonstrate  the 
impact  of  this  problem,  shadow  pixels  identified  by  the  automated  method  were 
compared  to  those  identified  by  hand  in  a  satellite  file  (see  Figure  8).  In  Figure  8,  panel 
A  represents  the  results  from  FASTEC,  where  green  identifies  shadow  pixels  and  yellow 
identifies  partial  shadow  pixels.  Due  to  resolution  differences  between  the  DSM  file  (one 
meter)  and  the  satellite  file  (2.4  meters),  the  resolution  from  the  automated  results  in 
Panel  A  was  matched  to  the  satellite  resolution  in  Panel  B.  This  led  to  between  four  and 
nine  DSM  pixels  comprising  each  pixel  in  Panel  A.  In  cases  where  every  DSM  pixel 
within  a  satellite  pixel  was  identified  as  shadow  by  FASTEC,  the  comparison  identified 
the  pixel  as  fully  in  shadow.  In  cases  with  a  mixture  of  shadow  and  sunlit  pixels  existed, 
the  comparison  identified  them  as  partial  shadow.  Panel  B  represents  the  results  of  a  by¬ 
hand  analysis  of  the  satellite  image.  Panel  C  shows  the  shadow  pixel  overlap  (162  pixels) 
between  Panels  A  and  B  without  perfonning  offset  correction.  By  shifting  the  latitude 
and  longitude  of  every  pixel  in  the  DSM  file  equally,  the  overlap  increased  to  218  pixels, 
as  shown  in  Panel  D.  Without  offset  correction,  38%  of  the  FASTEC  shadow  pixels 
were  actually  sunlit,  resulting  in  higher  shadow  radiance  and  higher  values  of  optical 
depth.  This  case  required  an  eight  meter  shift,  performed  by  hand,  to  maximize  the 
overlap.  Other  cases  required  shifts  of  100  meters  or  more,  most  likely  due  to  the 
imprecise  orthorectification  applied  to  the  satellite  files. 
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Figure  8.  Four-panel  plot  of  offset  between  satellite  and  DSM  files  for  Goddard  Space 
Flight  Center  on  October  26,  2006.  In  Panels  A  and  B,  green  is  shadow  and 
yellow  is  partial  shadow  pixels.  Panel  A  displays  the  results  from  FASTEC  and  B 
displays  the  results  from  a  by-hand  analysis  of  the  satellite  file.  Panels  C  displays 
overlap  between  A  and  B  without  offset  correction  (162  pixels).  Panel  D  displays 
overlap  with  offset  correction  (218  pixels). 


Regardless  of  the  amount  of  shift  required,  correcting  for  the  offset  between  the 
input  files  is  critical  to  FASTEC.  An  automated  offset  correction  procedure  that  adjusts 
the  latitude  and  longitude  of  the  upper  left  corner  in  the  DSM  file  to  find  the  lowest 
average  shadow  pixel  radiance  was  developed.  The  range  over  which  the  average 
radiances  are  computed  is  set  in  the  configuration  file.  Shadow  pixels  are  identified 
using  the  coherent  shadow  array,  with  correction  for  edge  contamination  (see  Chapter  III, 
Section  B.5.b.).  This  procedure  corrects  for  translation  error  between  the  DSM  and 
satellite  file.  It  does  not  correct  for  rotation  error  or  distortion  error.  To  improve 
computation  time,  the  correction  is  performed  with  a  coarse  and  fine  stage.  During 
coarse  correction,  the  coordinates  of  the  upper  left  corner  of  the  DSM  are  shifted  by  an 
increment  in  the  configuration  file,  nominally  4  pixels  in  both  the  north/south  and 
east/west  directions.  After  each  shift,  the  average  coherent  shadow  radiance  of  one 
satellite  channel  is  calculated  for  each  valid  target  shadow.  The  offset  that  provided  the 
lowest  average  radiance  is  the  starting  point  for  fine  correction.  During  fine  correction, 
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the  process  is  repeated  for  every  possible  shift  between  the  coarse  increment  and  the 
coarse  correction  best  fit.  After  the  lowest  average  shadow  radiance  is  found,  the  latitude 
and  longitude  of  the  DSM  file  is  adjusted  by  that  offset  interval. 

5.  Eliminate  Unwanted  Pixels 

Several  procedures  were  developed  to  eliminate  pixels  whose  radiances  are  not 
representative  of  the  radiance  in  shadow  or  sunlit  regions  because  of  reduced  diffuse 
radiance,  combined  shadow/sunlit  radiance,  and  differing  surface  reflectance.  Each 
procedure  is  geared  toward  testing  target  selection  theories  as  outlined  in  Vincent  (2006), 
Rivenbark  (2009),  and  Dombrock  (2007).  Figure  3  displayed  examples  of  output  from 
FASTEC  for  two  of  these  criteria.  Each  of  the  rejection  methods  are  described  in  more 
detail  in  the  remainder  of  this  section.  Different  variables  in  the  configuration  file  control 
execution  of  these  procedures,  thereby  making  them  easily  adjustable. 

a.  Pixels  That  Are  Too  Close  to  the  Building 

As  portions  of  the  sky  are  blocked  by  taller  structures,  the  amount  of 
diffuse  sky  radiance  reflecting  off  a  pixel  is  reduced,  thereby  lowering  radiance  values 
and  increasing  error  in  AOD  measurements.  Rivenbark  (2009)  and  Dombrock  (2007) 
suggested  selecting  shadow  pixels  far  from  the  shadow  generator  in  order  to  mitigate  the 
impact  of  this  effect.  FASTEC  includes  two  methods  of  eliminating  pixels  based  on  their 
proximity  to  the  shadow  generator;  the  first  eliminates  pixels  using  their  physical  distance 
from  the  shadow  generator  while  the  second  eliminates  a  percentage  of  the  pixels  within 
the  overall  shadow  (see  Figure  9). 
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Figure  9.  Identifying  pixels  by  proximity  to  shadow  generator.  hbidg  is  the  height  of  the 
building,  hb  is  the  height  of  the  line-of-sight  between  the  target  pixel  and  sun, 
distpix  is  the  distance  from  the  shadow  generator  to  the  target  pixel,  distmax  is  the 
farthest  extent  of  shadow,  and  puser  is  the  percentage  of  pixels  to  mark  invalid. 

In  the  interest  of  computational  efficiency,  these  calculations  are  actually 
performed  when  identifying  shadow  and  sunlit  pixels  (see  Chapter  III,  Section  B.2)  but 
are  stored  in  a  separate  array.  These  adjustments  are  separated  to  prevent  error  in  the 
offset  correction  algorithm.  Additionally,  it  allows  the  output  from  this  step  to  be  stored 
in  a  file  and  examined  separately  during  testing.  Assuming  the  radiance  is  lowest  near 
the  building  where  diffuse  sky  radiance  is  blocked  by  the  building,  performing  the  offset 
correction  after  this  step  would  choose  offset  values  that  place  the  shadow  region  too 
close  to  the  building,  the  exact  location  this  step  is  trying  to  avoid.  This  would  cause  the 
alignment  between  the  DSM  and  satellite  file  to  be  incorrect,  likely  resulting  in 
significantly  lower  sunlit  radiance  values  as  many  pixels  identified  as  sunlit  in  the  DSM 
would  fall  over  shadow  in  the  satellite  file. 

When  selecting  pixels  based  on  the  physical  distance  from  the  shadow 
generator,  the  distance  from  the  target  pixel  to  the  shadow  generator  is  compared  to  a 
threshold  as  defined  in  the  configuration  file.  If  the  distance  is  less  than  the  threshold,  the 
shadow  pixel  is  marked  invalid.  To  select  pixels  based  on  their  proximity  by  percentage 
of  the  overall  shadow,  the  ratio  displayed  in  Figure  9  is  compared  against  the  building 
percentage  rejection  threshold  in  the  configuration  file.  If  the  ratio  is  less  than  the 
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threshold,  the  pixel  is  marked  invalid.  The  results  from  this  step  are  integrated  into  the 
coherent  shadow  array  immediately  after  offset  correction  is  performed. 


b.  Pixels  Whose  Sky  Dome  Is  Blocked  by  Adjacent  Pixels 

As  described  in  the  preceding  section,  diffuse  sky  radiance  is  blocked  by 
any  tall  objects  near  a  pixel.  The  test  based  on  distance  from  the  building  addressed  the 
problem  between  shadow  pixels  and  the  shadow  generator.  To  fully  address  the  impact 
of  diffuse  sky  radiance,  the  automated  approach  must  examine  other  surrounding  objects. 
A  procedure  was  developed  to  examine  the  sky  dome  blockage  of  shadow  and  sunlit 
pixels  before  allowing  them  to  be  used  in  the  AOD  calculation.  The  sky  dome  is 
considered  the  region  of  the  sky  that  is  visible  to  the  target  pixel. 


Figure  10.  Celestial  dome  blockage.  The  dark  line  depicts  the  angle  representing  dome 

blockage  along  that  azimuth. 


The  procedure  determines  the  percentage  of  sky  dome  blockage  by 
averaging  the  blockage  along  equally-spaced  azimuths  (see  Figure  10).  The  angle 
between  every  pixel  within  the  search  range  and  the  target  pixel  is  computed  along  each 
azimuth  in  a  circle  around  the  target  pixel.  The  azimuth  interval,  distance  interval,  and 
search  range  are  assigned  using  the  configuration  file.  The  angle  ( 9  )  between  the  pixels 
is  computed  using  the  following  equation: 

f  \ 


9  =  tan  1 


h(r,,c2)-h(rtar,ctar) 

Jdy(r2-yj2+dx(c2-ctj\ 


180 

71 


(12) 
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where  (r2,C2)  and  (rtar,Ctar)  are  the  along-azimuth  pixel  and  target  pixel  respectively,  h(r,c) 
is  the  height  of  that  point,  and  dx  and  dy  are  the  pixel  resolution  in  the  north/south  and 
east/west  directions.  The  average  maximum  blockage  angle  along  each  azimuth  is 
divided  by  90  to  get  the  overall  blocked  percentage  of  the  sky  dome.  If  the  average 
blockage  exceeds  the  threshold  in  the  configuration  file,  the  pixel  in  the  coherent  shadow 
array  is  labeled  invalid.  The  output  from  this  procedure  is  integrated  into  the  coherent 
shadow  array  prior  to  offset  correction.  Additionally,  the  average  blockage  calculation 
results  are  stored  in  an  “average  blockage”  array. 

c.  Pixels  That  Are  Too  Close  to  the  Shadow  Edge 

Pixels  along  shadow  edges  are  partially  sunlit.  This  is  apparent  because 
the  radiance  values  increase  significantly  near  the  edge  of  the  shadow.  Figure  1 1 
displays  the  shadow  pixel  radiances  as  detennined  by  FASTEC.  In  this  example,  pixels 
near  the  center  of  the  shadow  have  spectral  radiance  values  between  15  and  17  RU  while 
radiances  get  larger  as  the  shadow  edge  is  neared.  This  pattern  is  reversed  for  sunlit 
pixels.  Including  these  edge  pixels  raises  shadow  radiance  and  lowers  sunlit  radiance, 
causing  computed  AOD  values  to  increase.  In  the  satellite  image,  one  would  expect  this 
contamination  to  affect  only  those  pixels  along  the  shadow  edge.  Due  to  the  nearest 
neighbor  method  of  retrieving  radiances  for  each  pixel  and  issues  with  offset  between  the 
images,  this  problem  may  impact  more  than  just  those  pixels  closest  to  the  shadow  edge. 
A  procedure  was  created  to  eliminate  contaminated  edge  pixels  from  the  computation. 
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66  68  69  63  45  21  18  17  18  18  18  18  18  24  31  26  35  48  49  45 

73  70  61  34  21  18  17  17  18  18  17  17  18  18  23  30  27  39  49  48 

68  69  67  54  23  19  18  18  17  17  18  17  17  17  17  17  17  23  31  28  34  47 

68  64  39  22  17  17  17  18  18  17|l6  17  17  17  16  17  17  18  21  27  25  38  47 

70  70  69  59  38  20  18  17  17  17  17  17|l7  16  16  16  16  16  16  17  17  17  20  31  25  37  46  47 

67  69  70  66  53  32  20  18  17  17  17  17  16  16  16  15  16  16  15  16  16  17  16  17  18  24  29  26  37  43  28  18 

69  69  64  54  26  19  18  18  17  17  17  16  16  16  16  15  16  16  16  16  16  16  16  16  17  18  25  28  24  26  19  17 

70  70  63  49  28  21  18  18  17  19  18  17  16  17  16  16  16  16  16  16  16  15  16  16  17  16  16  16|  17  25  29  19*17  16  16  16 

17  16  17  17  17  17  17  17  16  16  15  16  15  15  15  16  16  16  16  15  16  16  16  15  15  16  22  26  17  15  16  15 

18  16  16  16  17  16  16  16  16  15  16  15  15  16  16  16  15  16  16  15  15  15  15  15  15  15  16  18116  16  15  15  15 

55  50  27  17  16  16  16  16  16  15  15  15  15  15  15  15  15  16  16  15  15  15  15  15  15  15  15  15  15  14  15  15  15  15  15 

55  47  33  21  16  15  16  15  15  15  15  15  15  15  15  15  15  16  16  15  15  15  15  15  14  14  15  14  15  15  15  15  15  15  15  21  36 

57  44  36  37  20  16  16  16  16  15  15  15  15  15  15  15  15  15  15  15  15  15  15  15  15  14  15  15  15  15  16  15  15  15  16  20 

55  51  36  35  36  20  15  16  16  15  15  15  15  15  15  15  15  16  16  15  15  15  15  14  15  14  15  14  14  14  15  15  16  15  14  15  19  31 

57  48  34  39  31  18  15  15  15  16  15  15  15  15  15  15  15  15  15  15  15  15  15  15  15  15  14  15  14  15  15  15  15  14  15  15  23  32 

42  36  38  22  16  15  16  16  16  16  16  16  16  15  15  15  15  15  15  15  15  15  15  15  14  14  14  15  15  15  15  15  15  15  15  20 

45  30  18  17  16  17  17  17  16  16  16  15  15  15  15  15  15  15  14  15  15  14  14  14  15  15  15  16  15  14  15  15  15  27  37  / 

45  37  31  20  17  17  17  17  16  16  16  16  15  15  15  15  14  14  14  14  14  14  14  14  15  15  15  15  15  15  15  14  14  15 

38  35  33  23  18  17  17  16  16  16  15  15  15  15  15  14  14  15  14  14  14  14  14  14  14  15  15  15  14  14  14  15  15  18  26 

35  34  27  19  17  16  17  16  15  15  15  15  15  15  14  14  14  14  14  14  14  14  15  15  15  15  15  15  14  15  15 

35  34  30  22  17  16  15  15  15  15  15  14  15  15  14  14  15  14  14  14  14  14  15  15  14  14  15  14  15  15  17  31 

34  32  26  17  16  16  15  15  15  15  14  15  15  14  15  14  14  14  14  14  14  14  15  15  14  14  15  15  16  32  56 

34  34  30  18  16  15  15  15  15  15  14  15  14  14  14  14  14  14  14  15  14  14  15  14  14  15  16  17: 27  58  62 

34  34  32  22  16  15  15  15  15  15  15  14  14  14  14  14  14  14  14  14  14  14  14  14  14  16  17  22  54  61  55 

34  33  27  17  15  15  15  15  14  15  14  14  14  14  14  14  14  14  14  14  14  14  14  15  17  17  46  62  56 

34  34  30  18  15  15  15  15  15  15  14  14  14  15  14  14  14  14  14  14  14  16  18  32  62  4 

33  33  32  22  16  15  15  15  15  15  14  14  15  15  14  14  14  14  14  14  14  57  /  ^ 

33  32  26;15  15  15  15  15  14  15  15  14  14  14  14  14  14  15 

33  33  31  29  18  15  15  15  15  15  15  14  14  14  14  14  14  14  14  ^  * 

33  32  32  23  16  15  14  15  15  14  14  14  14  14  14  14  14 

35  34  34  27  17  15  15  15  14  14  14  14  14  14  14  14 

36  36  33  31  20  15  14  14  14  14  14  14  14  14  14 

37  33  32  24  16  15  14  14  14  14  14  14  14 

36  34  33  28  18  14  14  14  14  14  14  14 

32  33  33  31  21  15  14  14  14  14  14  14 

35  34  32  26IJ7  15  15  14  14  14  14 
35  35  33  32  21  15  15  14  14 
33  32  34  35  27®  15  15  14 
33  34  36  31  18ffl| 

35  36  37  35  / 


Figure  11.  Example  of  shadow  radiances  from  FASTEC.  Spectral  radiance  values  are 

e  2  11 

units  of  W  m"  sr‘  nm~  (RU).  Building  location  is  estimated. 


The  procedure  to  eliminate  shadow  edge  pixels  examines  every  pixel 


within  range  (as  assigned  by  the  configuration  file)  of  valid  shadow  pixels.  If  a  sunlit 
(shadow)  pixel  is  within  range  of  the  target  shadow  (sunlit)  pixel,  the  target  pixel  is 
labeled  invalid.  The  result  from  this  step  is  stored  in  an  “interior”  array  for  later  use. 
The  interior  array  is  integrated  into  the  coherent  shadow  array  prior  to  offset  correction  to 
reduce  the  shadows  used  in  offset  correction  to  those  most  likely  to  have  optical  depth 


calculations  performed. 


d.  Pixels  with  Elevation  Issues  or  Extreme  Radiances 


The  shadow  method  requires  the  shadow  and  sunlit  regions  to  have  the 
same  surface  reflectance.  Without  surface  reflectance  maps  of  comparable  resolution. 


this  requirement  is  impossible  to  fully  satisfy.  In  order  to  improve  the  likelihood  this 
requirement  is  met,  two  final  corrections  are  performed. 
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The  first  correction  examines  the  elevation  of  every  point  within  the 
shadow  region.  If  the  difference  between  the  maximum  and  minimum  elevation  exceeds 
a  threshold  in  the  configuration  file,  the  pixel  whose  elevation  falls  farthest  from  the 
mean  elevation  of  the  shadow  region  is  labeled  invalid.  This  process  is  repeated  until  the 
difference  between  the  max  and  min  elevation  no  longer  exceeds  the  threshold.  This 
procedure  eliminates  trees,  small  buildings,  cars,  or  other  tall  objects  that  most  likely 
have  different  surface  reflectance  within  the  shadow  region.  It  does  not  account  for 
different  surface  types,  such  as  asphalt  and  grass,  which  change  surface  reflectance.  A 
similar  procedure  for  sunlit  pixels  is  described  in  the  next  section. 

The  second  correction  is  the  only  correction  that  is  based  on  satellite  pixel 
radiance.  It  assigns  radiances  to  the  remaining  valid  shadow  pixels,  sorts  them  by 
radiance,  and  then  eliminates  percentages  of  the  extremes  based  on  thresholds  in  the 
configuration  file  (see  Figure  12).  The  pixels  are  sorted  by  radiance  as  they  are  read  into 
the  array  and  the  extremes  are  not  used  to  compute  the  average  radiance  of  the  shadow. 
This  procedure  eliminates  small  areas  in  the  shadow  region  with  different  surface 
reflectance  and  counteracts  some  offset  correction  error. 


Figure  12.  Histogram  of  radiance  values  with  extreme  radiance  rejection.  The  actual 
distribution  is  not  necessarily  normal  as  depicted  in  this  graphic. 
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6.  Identify  Accompanying  Sunlit  Pixels 


For  each  shadow  that  is  large  enough  (as  defined  in  the  configuration  file)  for 
AOD  computation,  sunlit  radiance  is  required.  The  procedure  to  compute  sunlit  radiance 
begins  by  associating  a  set  of  candidate  sunlit  pixels  based  on  their  proximity  to  the 
shadow.  The  search  region  for  valid  sunlit  pixels  is  limited  to  a  180  degree  region  around 
the  shadow,  perpendicular  to  the  sun  (see  Figure  13).  Using  an  azimuth  and  distance 
approach  similar  to  the  procedure  described  in  the  section  for  sky  dome  blockage,  this 
procedure  searches  for  candidate  sunlit  pixels. 


Non-shadow  Region 

Shadow 

pixels 


Figure  13.  Region  for  sunlit  pixel  search. 

After  identifying  the  initial  sunlit  candidates,  the  procedure  eliminates  unwanted 
pixels.  First,  sunlit  pixels  are  eliminated  based  on  their  sky  dome  blockage,  using  the 
same  procedure  described  in  Chapter  III,  Section  B.5.b.  Next,  it  eliminates  pixels  based 
on  elevation  using  a  similar  process  to  the  method  described  in  the  previous  section. 
Unlike  the  shadow  pixel  elevation  elimination,  this  process  compares  the  elevation  of 
each  sunlit  candidate  to  the  average  elevation  of  the  shadow  region.  Any  pixel  elevation 
that  falls  outside  the  tolerance,  as  defined  in  the  configuration  file,  is  eliminated  from  the 
sunlit  pixels  dataset.  Since  it  is  possible  to  associate  a  sunlit  pixel  with  more  than  one 
shadow,  this  elimination  is  done  for  each  individual  shadow.  This  approach  is  done  to 
align  the  shadow  and  sunlit  regions  in  an  effort  to  improve  confidence  that  they  have  the 
same  surface  reflectance. 
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The  final  elimination  step  removes  the  extremes  of  the  sorted  radiances  using  the 
same  approach  described  in  Chapter  III,  Section  B.5.d.  Finally,  the  average  sunlit 
radiance  is  computed  from  the  remaining  sunlit  pixels.  The  automated  shadow  detection 
program  described  above  greatly  increases  the  total  number  of  shadow/sunlit  pairs  for  a 
given  time  period  when  compared  to  manual  calculations.  The  last  major  section  of  the 
program  ties  in  past  research  to  this  automated  shadow  detection  algorithm. 

7.  Calculate  AOD  for  Each  Shadow 

The  AOD  calculation  is  performed  for  each  shadow/sunlit  pair  that  meet  the 
minimum  size  requirement  using  a  series  of  subroutines  to  compute  Equations  (1)  -  (4). 
The  subroutines  were  optimized  to  improve  computation  speed  but  follow  the 
methodology  used  in  past  research  (Vincent  2006,  Dombrock  2007,  Rivenbark  2009) 
with  corrections  to  the  issues  described  in  Chapter  II,  Section  B.4.  In  instances  where  the 
shadow  radiance  is  larger  than  the  sunlit  radiance  (due  to  poor  target  selection,  offset 
correction,  or  cloud  cover),  the  calculation  is  not  performed  and  dummy  TOD  and  AOD 
values  are  passed  to  the  output  file  to  allow  for  further  examination  of  those  shadows. 
The  output  from  the  automated  algorithm  is  described  in  detail  in  Chapter  III,  Section 
C.3. 

C.  DATA 


DSM-LIDAR  Satellite 


Figure  14.  Examples  of  input  data  used  in  this  research.  Both  images  are  over  the 
Goddard  Space  Flight  Center  (GSFC)  near  Washington  D.C.  The  QuickBird 
satellite  image  is  from  October  26,  2006. 
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1.  Digital  Surface  Model  (DSM) 


For  the  purposes  of  this  study,  the  term  Digital  Surface  Model  (DSM)  refers  to  the 
one  meter  resolution  elevation  maps  used  as  the  basis  of  the  automation  technique.  In 
these  maps,  each  data  point  (or  pixel)  contains  the  elevation  of  the  first  solid  object 
encountered  when  looking  down  at  the  surface  from  above.  Several  methods  have  been 
developed  for  creating  a  DSM.  Two  measures  are  commonly  employed  to  describe 
accuracy;  Circular  Error  of  90%  (CE90)  for  horizontal  accuracy  and  Linear  Error  of  90% 
(LE90)  for  vertical  accuracy.  CE90  is  the  minimum  diameter  of  the  horizontal  circle  that 
can  be  centered  on  all  photo-identifiable  Ground  Control  Points  (GCPs)  and  also  contain 
90%  of  their  respective  twin  counterparts  acquired  in  an  independent  geodetic  survey. 
LE90  is  the  linear  vertical  distance  that  90%  of  control  points  and  their  respective  twin 
matching  counterparts  acquired  in  an  independent  geodetic  survey  should  be  found  from 
each  other  (GeoVAR  2009).  Two  techniques  for  creating  high  resolution  DSM  data  are 
stereo-pairs  and  light  detection  and  ranging  (LIDAR). 

Stereo  pairs  consist  of  two  images  of  the  same  scene  taken  at  near  simultaneous 
times  with  slightly  different  view  angles.  Ground  control  points  provide  baseline 
location  and  elevation  infonnation  to  calibrate  the  image  and  improve  its  accuracy. 
DigitalGlobe  advertises  accuracy  of  6.5m  CE90  and  6.5m  LE90  (DigitalGlobe,  cited  June 
2010c)  although  Low  et  al.  (2006)  improved  CE90  to  5.0m  or  better. 

Another  method  of  creating  a  DSM  uses  LIDAR  instruments  flown  aboard 
aircraft.  The  U.S.  Army  Topographic  Engineering  Center  (TEC)  BuckEye  system 
employs  a  LIDAR  sensor  to  create  one  meter  resolution  DSM  data  of  0.5m  CE90  and 
0.3m  LE90  accuracy  (TEC  2009).  The  high  resolution  and  accuracy  of  this  data  makes  it 
a  superior  choice  for  the  automation  technique.  For  this  study,  the  LIDAR  data  is 
distributed  as  approximately  5km  x  5km  tiles  in  GeoTiff  file  format.  ENVI  4.5  is 
visualization  software,  capable  of  displaying  and  manipulating  data  in  a  multitude  of 
formats.  This  software  is  used  to  convert  these  DSM  files  into  flat  files  to  simplify  data 
ingest. 
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2. 


Satellite  Data 


High-resolution  commercial  satellite  data  from  the  IKONOS,  QuickBird  and 
WorldView-1  satellites  were  used  in  this  study.  QuickBird  offers  panchromatic  and 
multispectral  imagery  with  CE90  of  23m  and  horizontal  resolutions  of  60cm  and  2.4m, 
respectively.  WorldView-1  offers  a  panchromatic  imagery  with  CE90  of  6.5m  and 
horizontal  resolution  of  approximately  50cm  (DigitalGlobe,  cited  June  2010a/b). 
IKONOS  offers  multispectral  and  panchromatic  imagery  with  CE90  of  about  6m  and 
horizontal  resolutions  of  lm  and  3.28m,  respectively  (Dial  et  al.  2003,  IKONOS,  cited 
May  2010).  All  satellite  data  was  ordered  via  the  Internet  using  the  Web-based  access 
and  retrieval  portal  on  the  National  Geospatial-Intelligence  Agency  Web  site  in  NITF  2.0 
format.  This  study  used  basic  imagery  from  DigitalGlobe  (QuickBird,  WorldView-1). 
This  imagery  is  not  orthorectified  although  a  built-in  process  in  ENVI,  using  rational 
polynomial  coefficients,  was  used  in  an  effort  to  remove  some  distortion.  ENVI  software 
was  also  used  to  save  the  imagery  as  flat  files.  Although  Rivenbark  (2009),  found  that 
orthorectification  was  not  critical  when  implementing  the  shadow  technique  by-hand, 
FASTEC  needs  each  pixel  to  be  as  accurately  mapped  to  the  earth  as  possible,  and 
orthorectification  was  performed.  Using  fully  orthorectified  imagery  may  have 
significantly  improved  the  alignment  between  files  and  improved  the  results  (especially 
in  complex  urban  settings). 


21.6 


22.7 


17.4 


16.8 


Figure  15.  Nearest  neighbor  method.  In  this  example,  radiance  is  assigned  to  a  1  m 
resolution  DSM  pixel  (solid  box)  using  the  radiances  of  four  surrounding  2.4m 

satellite  pixels  (hatched  boxes). 
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Due  to  resolution  differences  between  the  DSM  and  satellite  file,  weighted 
averaging  was  used  to  assign  radiances  to  the  DSM  pixels  (see  Figure  15).  To  begin,  the 
latitude  ( cpD(r) )  and  longitude  (^D(r))  of  the  DSM  pixel  at  row  and  column  (r,c)  are 
determined: 

(PD(r)  =  (PDUL  -(r-l)dtpD  (13) 

^d(c)  =  ^Dul  +(c~l)dK  (14) 

where  the  coordinates  of  the  top  left  corner  ^cpDuL  ADul  )  and  change  in  latitude  and 

longitude  (dcpD,d)iD)  of  the  DSM  are  read  from  the  DSM  flat  file  header.  Next,  the 

latitude  and  longitude  are  converted  into  the  matching  row  and  column  of  the  satellite  file 
(rs,  cs)  using  the  following  equations: 


r  _(<PsUL-tPD(r))  ^ 

dtps 

c  _  (^p(c)~^sUL)  |  1 

dA,s 


(15) 

(16) 


where  the  coordinates  of  the  top  left  comer  ((pSuL  ,^SuL  j  and  change  in  latitude  and 


longitude  (dcps,dAs )  of  the  satellite  file  are  read  from  the  satellite  flat  file  header.  The 


satellite  row  and  column  number  are  real  numbers,  although  row  and  column  indexing 
uses  integers.  The  decimal  portion  of  these  numbers  is  used  to  determine  the  weight  of 
each  of  the  four  satellite  pixels  that  surround  the  DSM  pixel.  Following  Franke  and 
Nielson  (1980),  the  weights  are  assigned  using  inverse  weighting  based  on  that  decimal 
portion  using  variations  of  the  following  equations  (top  left  equations  shown): 


where  wtiL  is  the  weight  of  the  top  left  pixel,  R  is  a  constant  representing  the  maximum 
possible  fractional  distance  between  the  four  satellite  pixels.  Similar  equations  are  used 


33 


to  calculate  the  weights  of  the  top  right  (TR),  bottom  left  (BL),  and  bottom  right  (BR) 
pixels.  The  weighted  mean  of  the  four  satellite  pixels  relative  radiance  values  (DNtl, 
DNTr,  DNBl,  DNBr)  is  calculated  using: 

j-jN  _  (wtTLDNTL  +  wtTRDNTR  +  wtBLDNBL  +wtBRDNBR) 
wtTL+wtTR+wtBL  +  wtBR 

This  value  was  converted  to  absolute  spectral  radiance  (in  RU)  using  Equation  (8). 

Each  satellite  image  is  shipped  with  a  collection  of  files,  including  the  satellite  file 
and  a  metadata  file.  These  files  contain  critical  infonnation  about  the  location  of  the 
image,  the  channel  calibration  infonnation,  and  the  sun  and  satellite  geometry.  A 
procedure  was  written  to  read  the  key  data  from  these  metadata  files.  The  IKONOS 
metadata  files  are  fonnatted  differently  than  the  files  for  WorldView-1  and  QuickBird. 
They  also  lacked  the  necessary  calibration  infonnation.  To  conect  these  issues,  a 
program  was  written  to  reformat  the  IKONOS  metadata  file  into  a  suitable  fonn  and  to 
provide  channel-specific  calibration  information  that  was  obtained  from  GeoEye 
Corporation. 


3.  Output  Data 

FASTEC  generates  two  output  files  to  separate  general  infonnation  from  the 
AOD  calculations  and  simplify  statistical  analysis.  The  first  file,  identified  by  a 
“geninfo”  tag,  contains  general  information  about  that  particular  program  run.  This 
information  includes  the  satellite  name,  image  date  and  time,  the  satellite  resolution, 
number  of  channels,  the  sun  and  satellite  azimuth  and  elevation,  and  the  absolute 
calibration  factor,  effective  bandwidth,  solar  irradiance  and  5r  for  each  satellite  channel. 
The  “geninfo”  file  also  includes  the  name  and  resolution  of  the  DSM  file,  the  range  of  the 
DSM  over  which  the  program  was  run,  key  shadow/sunlit  adjustment  variables,  the 
number  of  shadows  identified,  the  DSM  offset  correction,  the  number  of  shadows  where 
AOD  calculations  were  performed,  and  a  line  stating  “Program  ended  successfully”  if  the 
program  completed  without  errors. 
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The  second  output  file,  identified  with  “AOD”  contains  the  AOD  calculation 
data  for  each  shadow  where  a  calculation  was  performed.  This  information  includes  a 
shadow  identification  number,  a  shadow  pixel  row  and  column  number  for  the  DSM  and 
satellite  file,  the  latitude  and  longitude  of  that  pixel,  and  pixel  counts  for  the  shadow  and 
sunlit  regions.  It  also  includes  the  channel  number,  average  shadow  and  sunlit  radiances, 
and  the  computed  TOD  and  AOD  value  for  each  satellite  channel. 

If  a  sub-region  is  less  than  or  equal  to  100  by  100  pixels,  the  user  may  create 
additional  output  files  by  setting  a  variable  in  the  configuration  file.  These  additional  text 
files  store  critical  array  data  used  in  the  automated  program.  They  are  identified  by 
“DSM”  and  a  number  as  described  in  Table  2.  If  sky  dome  blockage  is  not  tested,  file 
number  7  contains  zeroes.  File  numbers  8-11  are  updated  if  the  satellite  file  contains 
those  channels.  For  example,  if  the  satellite  file  contains  one  channel,  only  file  8  would 
be  updated  while  a  four  channel  file  would  update  files  8-11.  If  additional  satellite 
channels  were  available  (not  applicable  to  this  study),  additional  files  (numbered  beyond 
12)  would  be  generated  to  store  the  radiance  of  those  channels. 


ME  JIM 

Contents  of  file 

1 

Raw  DSM  elevation  data 

2 

Building  height 

3 

Distance  from  shadow  generator  with  satellite  geometry  correction 

4 

Same  as  3  with  pixels  that  are  too  close  shadow  generator  invalid 

5 

Same  as  3  with  shadow/ non-shadow  edge  pixels  invalid 

6 

Coherent  shadow  numbers  incorporating  invalid  pixels  files  4  and  5 

7 

Dome  block  %  for  each  valid  shadow/non-shadow  pixel 

8-11 

Absolute  spectral  radiance  for  each  channel  and  pixel 

Table  2.  List  of  contents  of  DSM  array  output  files. 


4.  Validation  Data 

Aerosol  Robotic  Network  (AERONET)  data  served  as  ground  truth  in  this 
research.  AERONET  is  a  worldwide  collection  of  robotic  sun  photometers  that  measure 
AOD  and  are  accurate  to  within  0.02  (Holben  et  al.  1998).  Three  levels  of  AERONET 
data  are  accessible  on-line  via  the  Goddard  Space  Flight  Center  Web  site.  Level  1.0  is 
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unscreened  data,  level  1.5  is  cloud-screened  data,  and  level  2.0  is  quality  assured  and 
cloud-screened.  Except  where  otherwise  noted,  level  2.0  data  was  used  for  validation. 

To  correct  for  wavelength  differences  between  AERONET  and  satellite  sensors, 
coefficients  for  the  third  order  polynomial  best-fit  curves  approximating  the  spectral 
distribution  were  determined  using  Microsoft  Excel.  The  coefficients  were  determined 
using  calculations  based  on  time-averaged  AOD  values  for  each  wavelength.  Although 
Excel  has  a  built-in  function  for  generating  polynomial  best  fit  curves,  an  error  in  the 
Excel  program  (detennined  by  comparing  the  built-in  plot  of  the  curve  with  a  manual  plot 
of  the  curve  using  the  Excel-identified  coefficients)  produces  the  wrong  coefficients  for 
polynomial  curves.  Therefore,  the  coefficients  were  determined  using  the  Excel  LINEST 
function  (Eng-tips,  cited  May  2010).  These  curves,  along  with  the  channel-specific 
spectral  response  curves  for  each  satellite,  were  used  to  calculate  sensor-specific  ground 
truth  AOD  using  Equation  (9).  IKONOS  spectral  response  data  is  available  on  the 
GeoEye  Web  site  while  QuickBird  and  WorldView-1  spectral  response  curves  are 
available  from  DigitalGlobe  upon  request.  A  third  order  polynomial  best  fit  curve 
estimated  AOD  values  to  within  0.003  of  truth  (see  Figure  16). 
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Figure  16.  Third  order  polynomial  best  fit  curve  with  raw  AERONET  data  points. 

Four  test  locations  were  chosen  based  on  the  availability  of  overlapping 
AERONET  stations  and  Buckeye  data.  The  four  stations  were  Tucson  (32N,  1 10W),  MD 


- Best  Fit  Curve 

.  AERONET 
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Science  Center  (39N,  76W),  University  of  Houston  (29N,  95W),  and  GSFC  (38N,  76W). 
Due  to  the  lack  of  high  AOD  cases  in  the  available  satellite  images  for  the  four  test 
locations,  an  additional  DSM  file  was  created  by  hand  near  the  Solar  Village  (24N,  46E) 
AERONET  station  (see  Figure  17).  The  heights  of  the  structures  were  estimated  by  using 
the  length  of  a  shadow  and  the  elevation  angle  of  the  sun  using  the  Solar  Village  image 
from  April  12,  2006  at  0754Z.  The  automated  program  was  designed  to  calculate  AOD 
from  the  two  satellite  dish  shadows.  The  building  was  included  to  ensure  the  automated 
method  did  not  include  those  pixels  as  part  of  the  sunlit  region.  The  satellite  dishes  are 
movable  objects,  and  orientation  changes  did  occur  between  images.  Orientation 
differences  are  noted  when  observed  in  testing. 


Figure  17.  Satellite  image  of  area  used  for  hand-generated  DSM  (image  captured  from 
Google  maps).  The  labeled  objects  were  the  only  structures  built  into  the  DSM. 

D.  COMPUTATIONAL  TECHNIQUES  AND  TESTING 

FASTEC  was  written  in  FORTRAN  with  OpenMP  for  parallel  processing.  Using 
parallel  processing  improved  the  computation  time  from  18.2  minutes  to  4.2  minutes  for 
the  entire  GFSC  DSM  file  using  two  Intel  Xeon  3.2  GHz  quad  core  CPUs  running 
LINUX  (see  Figure  18).  The  improvements  in  calculation  speed  achieved  through 
parallel  programming  made  it  possible  to  quickly  test  key  variables  across  ranges  of 
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possible  values.  Each  major  step  in  FASTEC  is  performed  using  a  different  subroutine  to 
allow  for  easy  modification  of  the  overall  process.  The  main  program  performs  a  few 
basic  tasks  and  calls  applicable  subroutines  as  needed. 
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Figure  18.  Plot  of  computation  time  vs.  number  of  processors.  Test  case  found  more 
than  50,000  shadows  in  25  km2  DSM  file  with  AOD  calculations  on  3015 

shadows. 


The  design  of  the  software  allowed  for  pixel-level  examination  of  the  automated 
approach  in  single  shadow  scenes.  Initial  testing  used  this  capability  to  study  the  impact 
of  different  shadow  and  sunlit  target  selection  criteria  and  led  to  baseline  values  for  AOD 
computations  with  every  available  image  and  DSM  file  to  test  the  overall  accuracy  of 
FASTEC.  A  separate  variable  combination  was  used  for  the  hand-generated  Solar 
Village  DSM  since  the  shadow  generators  were  comparatively  small.  For  example,  the 
satellite  dishes  in  Solar  Village  were  approximately  13  meters  tall  while  nine  shadow 
generators  in  GSFC  were  at  least  35  meters  tall. 
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IV.  RESULTS 


FASTEC  was  tested  to  ensure  that  it  properly  identified  building  shadows  and 
calculated  AOD.  Initial  testing  examined  the  primary  target  selection  criteria.  The 
results  of  these  tests  were  used  to  assign  configuration  values  to  perform  the  next  tests 
comparing  FASTEC  AOD  results  with  ground  truth  AERONET  data.  In  the  case  of 
multispectral  images,  different  channels  were  also  compared  with  each  other.  Finally,  the 
horizontal  resolution  of  the  DSM  was  reduced  to  examine  the  coarsest,  usable  resolution. 

A.  TARGET  SELECTION  CRITERIA  TESTS 

The  impact  of  key  target  selection  criteria  was  analyzed  using  three  satellite 
images;  a  QuickBird  image  from  October  26,  2006  at  1600Z  over  the  GSFC  DSM,  a 
WorldView-1  image  from  January  16,  2009  at  1600Z  over  the  GSFC  DSM,  and  a 
QuickBird  image  from  April  12,  2006  at  0754Z  over  the  hand-generated  Solar  Village 
DSM.  First,  pixels  along  the  shadow  edge  were  examined  to  identify  partial  shadow 
pixels.  After  establishing  an  edge  pixel  rejection  value  for  the  GSFC  building,  the 
minimum  shadow  generator  height  was  chosen  to  prevent  AOD  calculations  on  tree 
shadows.  This  was  followed  by  examining  the  effect  that  distance  from  the  shadow 
generator,  sky  dome  blockage,  and  extreme  radiance  rejection  had  on  computed  values  of 
AOD.  Except  where  otherwise  noted,  the  configuration  file  variables  were  assigned  the 
values  listed  in  Appendix  B  (Solar  Village)  and  C  (GSFC). 

1.  Rejection  by  Minimum  Shadow  Generator  Height  (GSFC) 

The  minimum  shadow  generator  height  was  detennined  for  the  GSFC  scene  after 
the  initial  single  shadow  edge  pixel  test  for  the  GSFC  building  (see  Chapter  IV,  Section 
A.2.b)  detennined  the  edge  pixel  rejection  variable  should  be  set  to  2.  For  this  test, 
FASTEC  was  run  using  the  October  26,  2006  GSFC  image  with  different  values  for  the 
minimum  shadow  generator  height.  Whenever  the  elevation  difference  between  a 
shadow  pixel  and  the  shadow  generator  was  below  this  threshold,  the  shadow  pixel  was 
marked  invalid.  Next,  each  shadow  in  the  AOD  output  file  was  located  in  the  satellite 
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image  by  hand  and  identified  as  being  caused  by  either  trees  or  buildings.  Figure  19 
displays  the  results  of  this  test  with  minimum  shadow  generator  height  values  between  25 
and  40.  Lower  values  were  not  analyzed  by  hand  due  to  the  large  number  of  shadows 
identified  and  the  fact  that  88.2%  of  the  shadows  were  tree  shadows  when  the  minimum 
shadow  generator  height  was  set  to  25.  With  settings  of  25  or  above,  the  number  of 
building  shadows  only  varied  by  two  buildings.  With  settings  of  35  or  40,  one  tree 
shadow  remained. 
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Figure  19.  Number  of  tree/building  shadows  vs.  minimum  shadow  generator  height  from 

October  26,  2006  GSFC  QuickBird  image. 

The  variation  in  AOD  values  was  also  examined  as  the  minimum  shadow 
generator  height  was  changed  (see  Figure  20).  In  this  graph,  the  50th  percentile  of  AOD 
was  plotted  for  different  settings  of  minimum  shadow  generator  height.  The  error  bars 
represent  the  25th  and  75th  percentile  of  AOD  values,  the  square  marks  the  50th 
percentile,  the  plus  symbols  mark  the  maximum  and  minimum  AOD  values,  and  the 
dotted  line  is  the  ground  truth  AERONET  AOD  value.  In  cases  where  the  minimum 
shadow  generator  height  is  less  than  30,  the  50th  percentile  value  and  range  is  larger  than 
cases  where  it  is  30  or  larger.  This  is  due  to  the  large  number  of  tree  shadows  in  the  low 
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minimum  shadow  generator  height  test  cases  because  tree  shadows  may  allow  some 
sunlight  to  pass  through  the  canopy,  increasing  the  shadow  radiance  and  the  computed 
AOD.  The  maximum  value  of  AOD  (not  plotted)  in  test  cases  below  30  falls  between 
0.96  and  2.754.  In  cases  where  the  number  of  building  shadows  exceeds  the  number  of 
tree  shadows,  the  range  of  AOD  values  between  the  25th  and  75th  percentile  shrinks 
from  approximately  0.10  to  less  than  0.02.  In  each  of  these  high  threshold  cases,  the 
maximum  is  due  to  the  one  tree  shadow  remaining  in  the  sample.  This  tree  shadow  can 
be  eliminated  by  increasing  the  minimum  pixel  threshold  to  70.  The  minimum  pixel 
threshold  sets  the  required  number  of  pixels  in  the  shadow  and  sunlit  region  for  AOD 
calculations  to  be  perfonned.  Removing  this  tree  shadow  from  the  sample  shrinks  the 
AOD  range  between  the  maximum  and  minimum  to  approximately  0.02,  similar  to  the 
range  between  the  25th  and  75th  percentile.  Due  to  these  findings,  tests  involving  the 
entire  GSFC  DSM  file  used  a  minimum  shadow  generator  height  of  35  and  single  shadow 
tests  on  the  GSFC  building  used  a  minimum  shadow  generator  height  of  15  to  capture  the 
approximately  20  meter  tall  GSFC  building.  However,  minimum  pixel  count  remained 
1 5  and  any  tree  shadows  were  eliminated  from  the  results  manually. 
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Figure  20.  Minimum  shadow  generator  height  vs.  50th  percentile  of  AOD.  Computed 
for  all  shadow  generators  from  the  October  26,  2006  GSFC  image.  Error  bars  are 
the  25th  and  75th  percentile,  plus  symbols  are  the  minimum  and  maximum  AOD 
(maximum  exceeds  0.35  in  some  cases),  and  the  dotted  line  is  ground  truth. 
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2.  Rejection  of  Edge  Pixels 


Hypothesis  I  states  that  it  is  important  to  avoid  partial  shadow  pixels  because  they 
increase  shadow  radiance  and  decrease  sunlit  radiance.  The  edge  pixel  rejection  variable 
eliminates  pixel  along  the  shadow/sunlit  edge.  Figure  21  shows  how  the  shadow  region 
for  the  north  satellite  dish  (left)  and  the  GSFC  building  (right)  change  when  the  edge 
pixel  rejection  variable  is  set  to  two.  In  this  figure,  yellow  represents  the  shadow 
generator,  dark  grey  represents  full  shadow  pixels  and  light  grey  represents  shadow 
pixels  that  were  eliminated  due  to  their  proximity  to  the  shadow  edge.  The  white  pixels 
between  the  shadow  and  shadow  generators  are  shadow  pixels  that  were  eliminated  due 
to  excessive  elevation  difference  or  because  they  were  blocked  from  view  of  the  satellite. 
For  example,  in  the  GSFC  scene,  pixels  on  the  northwest  side  of  the  shadow  region  were 
eliminated  because  they  were  blocked  from  view  of  the  satellite  while  some  of  the  pixels 
to  the  southeast  of  the  shadow  region  were  eliminated  due  to  elevation  difference  caused 
by  a  one  story  annex  off  that  wing  of  the  building. 


Figure  2 1 .  Example  of  edge  rejection  algorithm  on  Solar  Village  north  satellite  dish  (for 
April  12,  2006  image)  and  GSFC  building  (for  January  16,  2009  image).  Dark 
grey  pixels  are  full  shadow  pixels,  light  grey  are  rejected  due  to  proximity  to 
shadow  edge,  and  yellow  approximate  the  shadow  generator. 


a.  North  Satellite  Dish,  Solar  Village 

Figure  22  displays  the  average  pixel  radiance  based  on  their  proximity  to 


the  shadow  edge  for  the  north  satellite  dish  in  the  April  12,  2006  Solar  Village  image. 
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The  zero  line  marks  the  approximate  position  of  the  shadow  edge  so  that  every  pixel  falls 
on  one  side  of  the  shadow  edge  or  the  other.  The  average  radiances  of  the  pixels  closest 
to  the  shadow  edge  fall  between  more  interior  pixels.  Additionally,  the  standard 
deviation  of  pixels  that  are  fully  in  shadow  is  low  compared  to  the  standard  deviation  of 
the  edge  shadow  pixels.  This  is  partially  due  to  the  different  percentages  of  these  pixels 
in  direct  sunlight  but  position  error  between  pixels  in  the  DSM  and  satellite  image  also 
contributes.  Regardless,  including  these  edge  pixels  increases  the  range  of  radiances  in 
both  the  shadow  and  sunlit  region,  which  should  increase  the  uncertainty  of  the  computed 
AOD  from  one  shadow  to  the  next.  For  this  case,  setting  the  edge  pixel  rejection  variable 
to  one  eliminates  partial  shadow  pixels  without  reducing  the  number  of  usable  pixels 
below  the  arbitrarily  set  minimum  threshold  of  15.  Setting  the  interior  rejection  variable 
to  two  offers  the  lowest  uncertainty  of  the  shadow  region  but  it  decreases  the  number  of 
usable  shadow  pixels  62%. 
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Figure  22.  Average  pixel  radiance  versus  proximity  to  shadow  edge  for  the  north  satellite 
dish  with  the  April  12,  2006  Solar  Village  image.  The  x-axis  is  the  distance  from 
the  shadow  edge  (in  pixels)  with  zero  marking  the  shadow  edge,  error  bars  are 
one  standard  deviation,  and  plus  symbols  are  the  max  and  min  radiance. 
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When  calculating  optical  depth  for  this  shadow,  the  value  decreases  as  the 
edge  pixel  rejection  variable  increases  (see  Table  3).  This  decrease  is  most  significant  as 
the  edge  pixel  rejection  variable  changes  from  0  to  1.  In  this  example,  the  computed 
AOD  value  is  close  to  ground  truth  when  the  edge  pixel  rejection  variable  is  set  to  0 
because  the  partial  shadow  pixels  decrease  the  radiance  difference  by  over  1 1  RU  when 
compared  to  the  setting  of  1 .  As  before,  in  theory  these  pixels  have  some  direct  sunlight 
reaching  them  (partial  shadow)  and  the  apparent  improvement  in  accuracy  is  due  to 
misapplying  the  shadow  method.  For  the  remainder  of  tests  involving  the  Solar  Village 
scene,  the  edge  pixel  rejection  variable  was  set  to  1  to  minimize  the  uncertainty  of  the 
shadow  pixel  radiance  while  capturing  the  most  usable  shadows  across  all  available 
imagery.  The  extreme  pixel  rejection  variables  (Chapter  IV.A.5.b)  are  designed,  in  part, 
to  compensate  for  the  small  number  of  partial  shadow  pixels  that  remain  in  the  sample. 


Edge  Pixel 
Rejection 
Setting 

Shadow 

Pixel 

Count 

Non- 

Shadow 

Pixel  Count 

Shadow 

Radiance 

Non-shadow 

radiance 

TOD 

AOD 

Error 

0 

168 

1380 

80.98 

152.09 

0.303 

0.242 

-0.001 

1 

98 

1084 

75.44 

158.49 

0.255 

0.194 

-0.049 

2 

37 

775 

70.37 

160.89 

0.225 

0.164 

-0.079 

Table  3.  AOD  with  varying  edge  pixel  rejection  setting  for  the  north  satellite  dish 

in  the  April  12,  2006  Solar  Village  image.  Ground  truth  AOD  was  0.243. 


b.  Goddard  Space  Flight  Center  Building 

When  the  automated  method  was  applied  on  the  GSFC  building  pictured 
in  Figure  21,  errors  in  image  geo-referencing  and  solar  azimuth  and  elevation  cause 
minor  alignment  differences  between  the  DSM  and  satellite  file.  While  offset  correction 
improves  translation  error  between  the  images,  it  does  not  account  for  rotation  or 
distortion  errors.  In  addition,  the  presence  of  trees  to  the  southwest  of  the  building 
shadow  impacts  pixel  radiance  for  more  pixels  in  the  sunlit  region.  All  these  factors 
make  this  scene  more  complex  than  Solar  Village. 

Similar  to  the  Solar  Village  scene,  the  average  shadow  radiance  does  not 
stabilize  until  moving  three  pixels  in  from  the  shadow  edge  (see  Figure  23). 
Additionally,  sunlit  pixel  radiance  is  comparatively  low  when  within  two  pixels  of  the 
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shadow  edge.  Beyond  two  pixels  from  the  shadow  edge  the  average  radiance,  standard 
deviation,  and  difference  between  the  maximum  and  minimum  radiance  stabilize  in  the 
shadow  region.  Variation  in  the  sunlit  region  is  quite  high  due  to  the  significant 
contribution  of  the  surface  reflectance  dependent  direct  radiance  signal  to  the  total  signal. 
Based  on  this  examination,  it  appears  the  edge  pixel  rejection  variable  should  be  set  to  2. 
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Figure  23.  Average  pixel  radiance  vs.  proximity  to  shadow  edge  for  the  GSFC  building 
on  January  16,  2009.  The  x-axis  is  the  distance  from  the  shadow  edge  (in  pixels) 
with  zero  marking  the  shadow  edge,  error  bars  are  one  standard  deviation,  and 
plus  symbols  are  the  max  and  min  radiance. 

When  computing  AOD  while  varying  the  edge  pixel  rejection  variable  for 
the  same  scene,  fluctuations  are  quite  small  (see  Table  4).  Most  notable  is  the  change  in 
sunlit  radiance  when  transitioning  from  a  setting  of  one  to  two.  This  comparatively  large 
change  is  due  the  relatively  small  sunlit  pixel  count  and  the  tree  line  to  the  southwest. 
The  tree  line  shadow  causes  the  sunlit  region  to  shrink  on  several  sides  as  the  edge  pixel 
rejection  variable  is  increased,  eliminating  a  large  number  of  low  radiance  pixels.  This 
supports  setting  the  edge  pixel  rejection  variable  to  two  for  the  GSFC  scene. 
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Edge  Pixel 
Rejection 
Setting 

Shadow 

Pixel 

Count 

Non- 

Shadow 

Pixel  Count 

Shadow 

Radiance 

Non-shadow 

radiance 

TOD 

AOD 

Error 

0 

924 

357 

16.45 

34.43 

0.124 

0.054 

0.018 

1 

848 

256 

16.16 

35.75 

0.116 

0.046 

0.010 

2 

752 

180 

16.07 

36.77 

0.112 

0.042 

0.006 

3 

668 

137 

15.85 

36.75 

0.109 

0.039 

0.003 

4 

571 

101 

15.2 

36.81 

0.102 

0.032 

-0.004 

5 

460 

69 

14.8 

36.83 

0.098 

0.028 

-0.008 

6 

351 

40 

14.69 

36.57 

0.097 

0.027 

-0.009 

7 

250 

17 

14.6 

34.99 

0.100 

0.030 

-0.006 

Table  4.  AOD  with  varying  edge  pixel  rejection  setting  for  the  GSFC  building 

using  January  16,  2009  image.  Ground  truth  was  0.036. 


Next,  AOD  values  for  the  whole  GSFC  DSM  file  were  computed  for 
October  26,  2006  using  only  building  shadows  with  the  settings  described  in  Chapter  IV, 
Section  A.l.  The  October  26,  2006  image  was  chosen  over  the  January  16,  2009  image 
because  it  overlapped  the  entire  GSFC  DSM  region,  providing  nine  usable  shadows 
instead  of  the  four  usable  shadows  in  the  January  16  image.  For  these  computations, 
offset  correction  was  enabled,  the  minimum  shadow  generator  height  was  set  to  35,  the 
minimum  pixel  count  was  set  at  15,  and  tree  shadows  were  manually  eliminated  from  the 
results.  Figure  24  displays  the  results  as  the  edge  pixel  rejection  variable  is  changed  from 
zero  to  nine.  As  the  edge  pixel  rejection  variable  increased  from  0,  the  50th  percentile 
value  of  AOD  decreased.  Setting  the  edge  pixel  rejection  variable  to  two,  as  suggested 
by  the  previous  single  shadow  test,  selects  all  9  building  shadows  while  keeping  the 
range  of  computed  AOD  values  small.  Settings  beyond  5  eliminated  some  buildings 
from  the  sample. 
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Figure  24.  Aerosol  optical  depth  vs.  edge  pixel  rejection  setting  for  October  26,  2006 
GSFC  image.  The  black  box  is  the  50th  percentile,  the  error  bars  are  the  25th  and 
75th  percentile,  the  plus  symbols  are  the  max  and  min,  and  the  dashed  line  is 

ground  truth. 


3.  Rejection  Due  to  Distance  from  Shadow  Generator 

Hypothesis  II  states  that  it  is  important  to  choose  shadow  pixels  far  from  the 
shadow  generator  to  calculate  optical  depth.  Higher  elevation  objects  that  surround  a 
point  block  portions  of  the  sky,  reducing  diffuse  radiation  components  of  the  total 
radiance  signal  received  at  the  satellite  sensor.  Figure  25  shows  output  from  FASTEC 
when  a  percentage  of  pixels  near  the  generator  are  eliminated  from  the  total  shadow 
region  for  both  the  north  satellite  dish  in  Solar  Village  and  the  GSFC  building.  The  light 
grey  region  identifies  pixels  that  are  too  close  the  shadow  generator.  As  the  percentage 
increases,  more  of  the  total  shadow  region  is  kept.  Testing  of  the  importance  of  this 
parameter  was  done  using  these  scenes.  Additionally,  AOD  was  computed  for  the  GSFC 
DSM  file  as  the  building  shadow  percentage  threshold  was  varied  with  10%  increments. 
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Figure  25.  Example  of  target  pixel  rejection  based  on  proximity  to  the  shadow  generator. 
Yellow  indicates  building  pixels,  dark  grey  indicates  shadow  pixels,  and  light 
grey  indicates  eliminated  shadow  pixels.  Created  using  north  satellite  dish  in 
Solar  Village  on  April  12,  2006  (top  panels)  and  the  GSFC  building  on  January 

16,  2009  (bottom  panels). 

a.  North  Satellite  Dish,  Solar  Village 

Figure  26  depicts  the  spectral  radiance  based  on  distance  from  the  shadow 
generator  for  the  April  12,  2006  Solar  Village  image.  For  this  test  case,  the  edge  pixel 
rejection  variable  was  set  to  1  to  eliminate  most  of  the  partial  shadow  pixels  from  the 
sample.  In  general,  spectral  radiance  increases  as  the  distance  from  the  shadow  generator 
increases.  The  large  radiance  of  pixels  identified  as  being  within  1.3  meters  of  the 
satellite  dish  is  heavily  influenced  by  partial  shadow  pixels  in  the  southeast  comer  of  the 
shadow.  The  correlation  between  the  distance  from  shadow  generator  and  spectral 
radiance  is  about  0.45  when  including  every  shadow  pixel,  while  eliminating  the  closest 
pixels  from  the  sample  increases  correlation  to  0.77.  If  the  target  selection  variable  that 
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rejects  pixels  based  on  their  physical  distance  from  the  shadow  generator  were  set  to  2  or 
rejection  by  percentage  were  set  to  80%,  those  partial  shadow  pixels  would  be  marked 
invalid  and  eliminated  from  the  calculation. 
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Figure  26.  Spectral  radiance  vs.  distance  from  the  shadow  generator  (binned)  for  the 
Solar  Village  north  satellite  dish  on  April  12,  2006.  For  each  bin,  the  box  marks 
the  50th  percentile,  error  bars  are  the  25th  and  75th  percentile,  and  plus  signs  are 

the  maximum  and  minimum. 


Table  5  displays  the  results  for  AOD  calculations  on  the  north  satellite 
dish  when  varying  the  proximity  to  building  by  percentage  threshold  for  the  shadow 
region.  Generally,  AOD  increases  as  more  pixels  close  to  the  shadow  generator  are 
removed.  The  slight  AOD  decrease  between  70%  and  80%  is  due  to  higher  radiance 
from  partial  shadow  pixels  near  the  satellite  dish,  which  are  included  with  the  80% 
threshold.  The  average  shadow  radiance  increases  from  73.04  RU  to  78.06  RU  as  the 
percentage  of  building  shadow  is  decreased  from  80%  to  30%.  At  80%  or  greater,  the 
entire  shadow  is  used  in  the  computation.  Overall,  AOD  varies  by  0.025  across  the  range 
of  thresholds.  Fifty  percent  of  the  change  occurs  from  30  to  40%  when  the  pixel  count 


49 


nearly  doubles  and  the  average  shadow  radiance  decreases  from  78.06  to  75.43  RU.  In 


every  case,  the  computed  AOD  is  at  least  0.04  lower  than  ground  truth. 


%  of  shadow 

used 

Shadow 
pixel  count 

Non-shadow 
pixel  count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

100 

98 

1084 

73.04 

159.7 

0.179 

0.067 

90 

98 

1084 

73.04 

159.7 

0.179 

0.067 

80 

98 

1084 

73.04 

159.36 

0.180 

0.066 

70 

86 

1084 

72.89 

159.36 

0.179 

0.067  ! 

60 

73 

1084 

73.14 

159.7 

0.180 

0.066 

50 

55 

1063 

74.42 

159.48 

0.187 

0.059 

40 

38 

1021 

75.43 

160.09 

0.191 

0.055 

30 

21 

933 

78.06 

160.47 

0.204 

0.042 

Table  5.  Results  of  varying  proximity  to  building  threshold  percentage  for  north 

satellite  dish  from  April  12,  2006  Solar  Village  image.  Ground  truth  is  0.246. 


b.  Goddard  Space  Flight  Center  Building 

Figure  27  shows  the  radiance  for  pixels  binned  by  distance  from  the 
shadow  generator  for  the  GSFC  building  on  January  16,  2009.  In  general,  the  50th 
percentile  value  of  the  radiance  increases  slightly  with  distance  from  the  building.  The 
high  radiance  pixels  near  the  shadow  generator  appear  to  be  partial  shadow  pixels.  Part 
of  the  variation  in  spectral  radiance  across  the  shadow  is  due  to  variations  in  surface 
reflectance  from  trees,  cement,  and  grass  surfaces  in  the  scene.  Other  possible  reasons 
for  the  variations  across  the  shadow  include  rotation  or  distortion  error  between  the  DSM 
and  satellite  file  and  reflected  direct  sunlight  off  buildings  near  the  shadow.  The 
correlation  between  distance  from  the  shadow  generator  and  spectral  radiance  for  all 
shadow  pixels  is  0.20,  with  it  increasing  to  0.31  when  pixels  less  than  20  meters  from  the 
building  are  removed  from  the  sample.  This  is  much  lower  than  the  0.45  and  0.77 
correlations  for  the  Solar  Village  scene. 
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Distance  from  shadow  generator  (meters) 


Figure  27.  Spectral  radiance  vs.  distance  from  the  shadow  generator  (binned)  for  GSFC 
building  on  January  16,  2009.  For  each  bin,  the  box  marks  the  50th  percentile, 
error  bars  are  the  25th  and  75th  percentile,  and  plus  signs  are  the  maximum  and 

minimum. 


Figure  28  is  a  plot  of  pixel  radiances  in  the  shadow  region  from  the  same 
scene.  The  sun  is  located  to  the  SSE,  casting  the  shadow  to  the  NW.  Subjectively, 
radiance  appears  to  increase  with  distance  from  the  shadow  generator.  Additionally, 
radiance  is  highest  along  the  shadow  edge,  particularly  along  the  northern  edge  of  the 
shadow.  This  high  radiance  region  to  the  north  may  be  due  to  reflection  off  the  wing  of 
the  building  to  the  north  or  west.  In  both  this  scene  and  the  Solar  Village  scene  (not 
pictured),  there  is  a  region  in  the  center  of  the  shadow  that  appears  to  have  the  lowest 
radiance  values.  To  further  examine  this  relationship,  another  shadow  in  the  GSFC  DSM 
(pictured  in  Figure  25)  was  examined.  The  shadow  cast  by  this  tall,  isolated  building  was 
much  longer  than  wide.  In  this  case,  the  correlation  between  distance  from  the  shadow 
generator  and  spectral  radiance  was  only  0.09  (results  not  pictured). 
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Figure  28.  Plot  of  spectral  radiance  by  pixel  for  the  GSFC  building  from  January  16, 
2009.  The  GSFC  building  surrounds  all  sides  except  the  southwest.  The  sun  is 

located  to  the  southeast. 

Table  6  displays  computed  AOD  values  for  the  GSFC  building  as  the 
building  shadow  percentage  threshold  changes  from  100  to  20%  in  10%  increments.  For 
this  test,  offset  correction  was  enabled  and  the  edge  pixel  rejection  threshold  was  set  to 
two.  Generally,  AOD  increased  as  the  threshold  decreased  and  limited  the  shadow  region 
to  pixels  further  from  the  building.  The  overall  variation  in  AOD  is  lower  than  the  Solar 
Village  case  partly  because  of  the  larger  number  of  shadow  pixels  associated  with  this 
building  (667)  and  the  large  width  of  the  shadow  generator  (about  23  pixels)  relative  to 
the  length  (about  15  pixels).  Otherwise,  the  GSFC  building  and  Solar  Village  scene  show 
increasing  AOD  due  to  increased  shadow  radiance  as  more  of  the  pixels  near  the  building 
are  rejected. 
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%  of  shadow 

used 

Shadow 
pixel  count 

Non-shadow 
pixel  count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

100 

751 

179 

15.65 

37.92 

0.034 

0.002  | 

90 

751 

179 

15.65 

37.92 

0.034 

0.002 

80 

733 

179 

15.64 

37.92 

0.034 

0.002 

70 

626 

179 

15.74 

37.92 

0.035 

0.001 

60 

501 

167 

15.9 

37.82 

0.037 

-0.001 

50 

388 

143 

16.11 

37.52 

0.040 

-0.004  | 

40 

276 

101 

16.41 

37.62 

0.043 

-0.007 

30 

168 

71 

16.77 

38.84 

0.044 

-0.008 

20 

68 

28 

17.28 

44.03 

0.035 

0.001 

Table  6.  Results  of  varying  proximity  to  building  threshold  percentage  for  GSFC 

building  wing  from  January  16,  2009  GSFC  image.  Ground  truth  is  0.036. 


Figure  29  displays  the  range  of  AOD  values  for  the  entire  GSFC  DSM 
using  the  October  26,  2006  image.  For  this  test,  offset  correction  was  enabled  and  the 
building  shadow  percentage  was  varied  from  10%  to  100%  with  a  10%  increment.  For 
thresholds  above  50%,  a  lone  tree  shadow  with  high  AOD  (greater  than  0.25)  was 
eliminated  manually.  For  thresholds  50%  and  below,  tree  shadows  were  not  present. 
Otherwise,  this  graph  shows  minor  variation  in  both  the  50th  percentile  and  range  of 
AOD  values  for  all  shadows  in  the  scene.  In  every  case,  the  result  shows  a  low  bias  of 
approximately  0.05  and  the  spread  across  all  AOD  values  varied  between  0.026  and 
0.055.  These  results  do  not  show  the  percent  of  building  shadow  threshold  as  particularly 
important  to  the  shadow  method  although  the  low  optical  depth  of  the  scene  may 
contribute  to  the  lack  of  change  in  AOD.  For  scenes  other  than  Solar  Village,  the 
building  shadow  percentage  threshold  will  be  set  to  50%  to  increase  the  average  radiance 
slightly.  To  keep  the  Solar  Village  shadow  pixel  counts  above  the  minimum  pixel 
threshold  and  reduce  partial  shadow  pixels  near  the  shadow  generator,  the  physical 
distance  from  the  shadow  generator  threshold  will  be  set  to  two. 
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Figure  29.  AOD  vs.  percent  of  building  shadow  used  for  the  October  26,  2006  GSFC 
image.  The  black  box  is  the  50th  percentile,  the  error  bars  are  the  25th  and  75th 
percentile,  the  plus  symbols  are  the  max  and  min,  and  the  dashed  line  is  ground 

truth. 


4.  Rejection  Due  to  Sky  Dome  Blockage 

Hypothesis  III  states  that  sky  dome  blockage  is  important  when  choosing  valid 
shadow  and  sunlit  pixels  to  calculate  AOD.  This  is  due  to  blocked  diffuse  sky  radiance 
(negative  contribution  to  pixel  radiance)  and  direct  reflection  off  taller  objects  (positive 
contribution  to  pixel  radiance).  The  sky  dome  blockage  threshold  was  varied  on  the 
January  16,  2009  image  over  the  GSFC  building  to  verify  the  sky  dome  blockage 
procedure  worked  properly.  Figure  30  displays  the  change  to  the  shadow  region  with  the 
sky  dome  blockage  threshold  set  to  30%.  In  this  example,  light  grey  signifies  invalid 
shadow  pixels,  which  exceeded  the  sky  dome  blockage  threshold,  dark  grey  identifies 
valid  shadow  pixels,  and  yellow  approximates  the  GSFC  building.  The  invalid  region  is 
due  to  higher  sky  dome  blockage  from  the  building  wing  to  the  northeast  of  the  shadow 
and  to  a  lesser  extent  the  wing  to  the  northwest. 
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Figure  30.  Example  of  sky  dome  blockage  for  the  January  16,  2009  image  of  the  GSFC 
building.  Yellow  is  building  pixels,  dark  grey  is  valid  shadow  pixels,  and  light 
grey  is  invalid  shadow  pixels  from  30%  or  greater  sky  dome  blockage. 


a.  North  Satellite  Dish,  Solar  Village 

The  north  satellite  dish  in  Solar  Village  provided  a  simple  scene  to  test  the 
significance  of  sky  dome  blockage  as  it  relates  to  spectral  radiance.  The  presence  of  the 
small  building  to  the  west  of  the  dish  impacts  sky  dome  blockage  and  provides  a  way  to 
test  the  accuracy  of  the  calculation.  The  maximum  vertical  blockage  of  the  4  m  tall 
building  is  about  42  degrees  for  the  closest  shadow  pixels.  The  building  affects  up  to  13 
of  the  36  azimuths  used  in  the  computation  with  the  azimuth  interval  set  to  10  degrees. 
Therefore,  the  maximum  contribution  of  this  building  to  sky  dome  blockage  within  the 
shadow  is  about  0.119  (see  Figure  31).  By  subtracting  the  sky  dome  blockage  values 
with  the  building  removed  from  the  blockage  with  the  building  included,  the  maximum 
difference  was  0.1 19  and  decreased  as  distance  from  the  building  increased.  Considering 
the  typical  sky  blockage  of  shadow  pixels  within  this  shadow  is  about  0.30,  the  maximum 
contribution  of  the  building  to  sky  dome  blockage  is  40%.  The  correlation  between  sky 
dome  blockage  and  distance  from  the  shadow  generator  is  -0.79  with  the  building  and  - 
0.87  without  the  building.  Correlation  is  not  -1.0  because  of  the  presence  of  the  satellite 
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dish  to  the  south,  variation  in  the  north  satellite  dish  pixel  elevations,  and  possible 
incomplete  sampling  angles  in  the  sky  dome  blockage  calculation. 


Figure  3 1 .  Sky  dome  blockage  test  using  April  12,  2006  Solar  Village  image.  Grey 
region  is  shadow  pixels  identified  by  the  automated  method.  Numbers  mark 
geometrically  relevant  information  for  the  calculation  on  the  highlighted  pixel. 

Figure  32  is  a  graph  of  the  spectral  radiance  versus  sky  dome  blockage,  in 
2.5%  bins,  for  shadow  and  sunlit  pixels  for  the  north  satellite  dish.  Partial  shadow  pixels, 
because  of  the  edge  pixel  rejection  variable  being  one,  explain  the  large  maximum 
radiance  values  in  some  of  the  shadow  bins.  Regarding  shadow  pixels,  radiance 
decreases  with  increasing  dome  blockage  until  the  40%  bin  where  the  radiance  steadily 
increases.  The  increase  may  be  due  to  diffuse  radiance  from  the  4  m  tall  building,  partial 
shadow  pixels  along  the  shadow/shadow  generator  edge,  or  reflection  of  direct  sunlight 
off  the  side  of  the  building  into  the  shadow.  When  the  sky  dome  blockage  threshold  for 
the  shadow  region  was  varied  from  30%  to  50%,  AOD  changed  by  less  than  0.01  (graph 
not  pictured).  Dropping  the  threshold  below  30%  provided  less  than  the  15  valid  shadow 
pixels  necessary  to  perform  the  calculation. 
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The  variation  in  the  radiance  of  the  sunlit  pixels  in  Figure  32  appears 
random  with  no  discernable  trends  as  sky  dome  blockage  increases.  This  is  likely  due  to 
the  overwhelming  influence  of  the  direct  portion  of  the  total  signal  radiance.  The 
relatively  low  minima  between  15  and  25%  dome  blockage  were  due  to  five  partial 
shadow  pixels  along  the  edge  of  the  building  to  the  west  of  the  satellite  dish.  AOD  varied 
by  less  than  0.003  as  the  sunlit  pixel  sky  dome  blockage  maximum  threshold  increased 
from  5%  to  100%.  Based  on  these  results,  sky  dome  blockage  does  not  appear  to  be  very 

important  for  either  the  shadow  or  sunlit  region  in  this  Solar  Village  image. 
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Figure  32.  Spectral  radiance  of  sky  dome  blockage  bins  for  the  April  12,  2006  Solar 
Village  image.  For  each  bin,  the  box  (shadow)  or  diamond  (sunlit)  is  the  50th 
percentile,  error  bars  are  the  25th  and  75th  percentile,  plus  symbols  are  the 
maximum,  and  “x”  symbols  are  minimum  radiances. 
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b.  Goddard  Space  Flight  Center  Building 


The  Goddard  Space  Flight  Center  building  scene  provided  real  DSM  data 
to  analyze  the  impact  of  sky  dome  blockage  on  shadow  and  sunlit  regions.  The  20  meter 
tall  building  is  composed  of  several  wings  off  a  curved  connecting  building  (see  Figure 
33).  The  shadow  region  is  northwest  of  the  southernmost  wing  with  the  sunlit  region 
west  of  the  shadow.  The  least-blocked  shadow  pixels  are  in  the  southwest  portion  of  the 
shadow  while  the  most-blocked  shadow  pixels  are  along  the  northeast  edge.  The  least- 
blocked  sunlit  pixels  are  on  the  southeast  side  while  the  most-blocked  shadow  pixels  are 
along  the  north  side  of  the  sunlit  region. 


Figure  33.  Shadow  used  for  sky  dome  blockage  testing  for  Goddard  Space  Flight  Center 

building,  January  16,  2009  image. 

Figure  34  is  a  graph  of  the  spectral  radiance  as  sky  dome  blockage,  in  4% 
bin  intervals,  varies  for  shadow  and  sunlit  pixels  for  the  GSFC  building  using  the  January 
16,  2009  image.  For  the  shadow  bins,  the  maximum  radiance  pixels  fall  along  the  outer 
edge  of  the  shadow  suggesting  they  may  be  partial  shadow  pixels  caused  by  position 
error,  despite  offset  correction  being  performed,  between  the  DSM  and  satellite  file.  As 
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sky  dome  blockage  percentage  increases,  radiance  decreases  slightly,  only  to  increase 
above  32%.  For  shadow  pixels,  correlation  between  sky  dome  blockage  percentage  and 
spectral  radiance  is  0.35. 

As  blockage  increases  for  sunlit  pixels,  the  spectral  radiance  also  increases 

but  144  of  the  179  sunlit  pixels  are  less  than  16%  blocked.  The  lack  of  blocked  pixels 

above  16%  helps  explain  the  large  changes  in  the  50th  percentile  value  of  spectral 

radiance  at  those  blocked  percentages.  Correlation  between  sky  dome  blockage  and 

spectral  radiance  is  0.41  in  the  sunlit  region.  As  discussed  previously,  when  examining 

sunlit  pixels,  the  direct  reflection  portion  of  the  total  signal  dominates  and  it  is  difficult  to 

determine  the  impact  of  sky  dome  blockage.  If  it  were  possible  to  remove  the  direct 

reflectance  portion  of  the  signal,  a  pattern  might  emerge.  Otherwise,  radiances  in  the 

sunlit  region  are  highest  near  the  building.  This  is  likely  due  to  higher  surface  reflectance 

for  those  pixels  although  additional  radiance  from  reflection  of  direct  sunlight  off  the 

building  may  contribute. 
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Figure  34.  Spectral  radiance  of  sky  dome  blockage  bins  for  the  GSFC  building  on 
January  16,  2009.  For  each  bin,  the  box  (shadow)  or  diamond  (sunlit)  is  the  50th 
percentile,  error  bars  are  the  25th  and  75th  percentile,  plus  symbols  are  the 
maximum,  and  minus  symbols  are  the  minimum  radiances. 
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Table  7  displays  the  results  of  FASTEC  as  the  shadow  sky  dome  blockage 
maximum  threshold  is  varied.  Although  the  AOD  values  are  quite  small  in  each  case,  the 
average  shadow  radiance  and  AOD  increase  as  the  sky  dome  blockage  threshold 
increases.  Part  of  the  increase  is  due  to  the  higher  radiances  associated  with  high  sky 
dome  blockage  pixels  although  this  increase  is  partially  countered  by  variation  in  offset 
correction  as  the  threshold  changes  (see  Table  9).  The  change  in  offset  correction  values 
also  contributes  to  the  variation  in  sunlit  radiance.  If  offset  correction  perfectly  aligned 
the  DSM  and  satellite  files,  sunlit  radiance  would  vary  much  less  than  it  does  in  this 
table. 


snaaow  sKy 
dome  blockage 
threshold 

Shadow 
pixel  count 

Non-shadow 
pixel  count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

0.15 

20 

152 

14.09 

34.31 

0.026 

0.20 

184 

179 

14.42 

34.66 

0.028 

-0.008 

0.25 

324 

180 

14.75 

37.6 

0.025 

-0.011 

0.30 

411 

180 

14.81 

37.6 

0.026 

-0.010 

0.35 

509 

180 

15.03 

38.39 

0.026 

-0.009 

0.40 

614 

179 

15.06 

37.72 

0.028 

-0.008 

0.45 

708 

179 

15.36 

37.92 

0.031 

-0.005 

0.50 

750 

179 

15.65 

37.92 

0.034 

-0.002 

0.55 

751 

179 

15.65 

37.92 

0.034 

-0.002 

Table  7.  Results  of  varying  shadow  sky  dome  blockage  maximum  threshold  for 

January  16,  2009  image  of  GSFC  building  wing.  Ground  truth  is  0.036. 


Table  8  shows  the  results  from  FASTEC  as  the  sunlit  sky  dome  blockage 
maximum  threshold  is  varied  for  the  GSFC  building  shadow  using  the  January  16,  2009 
image.  Generally,  sunlit  radiance  increases  as  the  sunlit  sky  dome  blockage  threshold  is 
increased  as  expected  from  Figure  34.  Since  shadow  radiance  is  fixed,  AOD  decreases 
but  the  overall  change  is  quite  small.  This  is  due  to  the  majority  of  sunlit  pixels  having 
sky  dome  blockage  below  15%. 
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Non-snaaow  sKy 
dome  blockage 
threshold 

Shadow 
pixel  count 

Non-shadow 
pixel  count 

Shadow 

Radiance 

Non-shadow 

Radiance 

AOD 

Error 

0.05 

751 

22 

15.65 

36.87 

0.037 

0.001 

0.10 

751 

96 

15.65 

37.37 

0.035 

-0.001 

0.15 

751 

137 

15.65 

36.66 

0.037 

0.001 

0.20 

751 

156 

15.65 

36.93 

0.036 

0.001 

0.25 

751 

169 

15.65 

37.42 

0.035 

-0.001 

0.30 

751 

173 

15.65 

37.52 

0.035 

-0.001 

0.35 

751 

178 

15.65 

37.83 

0.034 

-0.002 

0.40 

751 

179 

15.65 

37.92 

0.034 

-0.002 

0.45 

751 

179 

15.65 

37.92 

0.034 

-0.002 

Table  8.  Results  of  varying  the  sunlit  sky  dome  blockage  maximum  threshold  for 

January  16,  2009  image  of  GSFC  building  wing.  Ground  truth  is  0.036. 


Finally,  the  sky  dome  blockage  thresholds  were  varied  for  the  entire  GSFC 
DSM  file  using  the  October  26,  2006  image.  Figure  35  is  a  plot  of  the  results  of  running 
FASTEC  while  varying  the  sunlit  sky  dome  blockage  maximum  threshold.  This  plot 
shows  little  variation  in  AOD  calculations  as  the  threshold  changes.  This  is  due  to  the 
fact  that  almost  all  the  pixels  in  the  sunlit  regions  have  sky  dome  blockage  less  than  0.05 
as  evidenced  by  the  average  sunlit  pixel  count  increasing  from  563  to  677  as  the 
threshold  increases  from  0.05  to  0.30.  At  thresholds  beyond  0.30,  sunlit  pixel  count  does 
not  change. 
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Figure  35.  AOD  vs.  sunlit  pixel  sky  dome  blockage  maximum  thresholds  for  the  GSFC 

DSM  file  from  January  16,  2009.  For  each  threshold,  the  box  is  the  50th 
percentile,  the  error  bars  are  the  25th  and  75th  percentile,  and  plus  symbols  mark 

the  maximum  and  minimum  AOD. 

Figure  36  is  a  plot  of  the  results  of  running  FASTEC  while  varying  the 
shadow  sky  dome  blockage  maximum  threshold.  The  minimum  AOD  (-0.11)  is  below 
the  y-axis  for  the  0.05  threshold.  Setting  shadow  sky  dome  blockage  thresholds  to  more 
than  0.60  does  not  eliminate  any  pixels.  As  opposed  to  changing  sunlit  dome  blockage 
maximum  threshold,  despite  small  AOD  values,  changes  in  the  50th  percentile  values  are 
apparent.  These  changes  mirror  changes  in  the  row  and  column  offset  for  the  different 
thresholds  (see  Table  9).  Where  the  row  and  column  offset  values  are  similar,  the  results 
in  the  figure  are  similar.  For  example,  when  the  threshold  is  30%  or  higher,  offset 
correction  adjustments  are  within  two  pixels  of  the  baseline  and  the  difference  between 
the  maximum  and  minimum  AOD  is  less  than  0.002.  Based  on  the  inconclusive  results 
of  these  tests,  sky  dome  blockage  thresholds  will  be  disabled  for  later  testing. 
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Figure  36.  AOD  vs.  shadow  pixel  sky  dome  blockage  maximum  thresholds  for  the  GSFC 
DSM  file  from  January  16,  2009.  Boxes  are  the  50th  percentile,  the  error  bars  are 
the  25th  and  75th  percentile,  and  +  symbols  mark  the  max  and  min  AOD. 


Shadow  sky 
dome  blockage 
threshold 

Row 

Offset 

Col 

Offset 

Distance  from 
baseline  (2,9) 

0.05 

13 

11 

11.18 

0.10 

10 

9 

8.00  [ 

0.15 

10 

9 

8.00 

0.20 

9 

9 

7.00 

0.25 

6 

8 

4.12 

0.30 

4 

9 

2.00 

0.35 

3 

8 

1.41  i 

0.40 

3 

8 

1.41 

0.45 

3 

9 

1.00  j 

0.50 

1 

8 

1.41 

0.55 

2 

9 

0.00  | 

Table  9.  Offset  correction  values  for  different  shadow  sky  dome  blockage 

thresholds  for  January  16,  2009  image  over  the  GSFC  DSM.  Baseline  offset 
correction  determined  by  disabling  the  sky  dome  blockage  maximum  threshold. 
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5. 


Rejection  of  Extreme  Radiances 


a.  North  Satellite  Dish,  Solar  Village 

Hypothesis  IV  states  that  eliminating  the  extreme  radiances  values  from 
the  sample  improves  AOD  calculation.  Figure  37  shows  the  sorted  spectral  radiances  of 
shadow  and  sunlit  pixels  for  the  north  satellite  dish  in  Solar  Village  from  the  April  12, 
2006  image  using  FASTEC.  The  first  approximately  100  sunlit  pixels  increase  in 
radiance  steeply  due  to  partial  shadow  pixels  along  the  northwest  side  of  the  small 
building  west  of  the  satellite  dish  and  partial  shadow  pixels  from  both  satellite  dish 
shadows.  Including  these  pixels  in  the  sunlit  region  average  causes  a  slight  decrease  in 
the  average  radiance  for  the  region  although  the  large  number  of  sunlit  pixels  lessens  the 
influence  of  these  darker  pixels.  When  the  105  darkest  (lowest  10%)  pixels  are  removed 
from  the  sample,  the  average  sunlit  radiance  increases  by  1.4  RU,  however,  when  the 
darkest  and  brightest  extremes  are  removed  equally  across  the  sunlit  region,  the  change  in 
radiance  is  quite  small  (see  Table  11). 

For  shadow  pixels,  the  overall  slope  of  the  increasing  radiances  is  steeper 
and  the  six  brightest  pixels  have  noticeably  higher  radiances  than  the  rest  of  the  sample. 
The  overall  steeper  slope  of  the  shadow  region  when  compared  with  the  sunlit  region 
appears  to  be  partly  due  to  the  smaller  number  of  pixels  in  the  shadow  region  (98  vs. 
1059  pixels).  The  significantly  brighter  pixels  are  partial  shadow  pixels  along  the 
shadow  edges.  Removing  these  six  pixels  causes  average  radiance  to  decrease  by  almost 
1.2  RU.  Removing  these  extreme  radiance  pixels  from  both  regions  results  in  a  2.6  RU 
increase  in  the  radiance  difference  and  computed  AOD  to  be  .0083  lower.  Although  this 
lowering  of  AOD  pushes  the  results  further  from  ground  truth,  it  eliminates  invalid  sunlit 
pixels  and  properly  applies  the  shadow  method. 

When  the  darkest  and  brightest  extremes  are  removed  equally  from  the 
shadow  region,  the  largest  change  in  shadow  radiance  occurs  when  10%  of  the  extremes 
are  removed  (see  Table  10)  and  radiance  difference  increases  from  86.8  RU  to  87.6  RU, 
decreasing  AOD  from  0.192  to  0.171.  The  radiance  difference  change  is  less  than  0.2 
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RU  as  additional  percentages  of  the  extremes  are  eliminated.  For  sunlit  pixels,  radiance 
difference  changes  by  0.3  or  less  across  all  combinations  (see  Table  1 1). 


0  100  200  300  400  500  600  700  800  900  1000  1100 


Sorted  pixels,  Solar  Village  04122006 


Figure  37.  Spectral  radiance  of  usable  pixels  for  shadow  and  sunlit  regions  for  Solar 
Village  north  satellite  dish  and  the  April  12,  2006  satellite  image. 


Min  pixel 
position 

Max  pixel 
position 

Min  pixel 
radiance 

Max  pixel 
radiance 

1 

98 

66.3 

96.8 

9 

88 

67.9 

77.9 

19 

78 

68.9 

75.7 

29 

68 

70.3 

74.0 

39 

58 

71.1 

73.4 

49 

49 

72.1 

72.1 

Radiance 

difference 


86.8 


87.6 


87.6 


87.7 


87.7 


87.8 


Table  10.  Results  of  varying  extreme  radiance  thresholds  for  Solar  Village  north 
satellite  dish  shadow  region  and  the  April  12,  2006  satellite  image. 
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19 

mi 

Min  pixel 
position 

Max  pixel 
position 

Min  pixel 
radiance 

Max  pixel 
radiance 

Average  non¬ 

shadow 

radiance 

Radiance 

difference 

0 

100 

1 

1059 

115.3 

179.3 

159.9 

86.8 

10 

90 

105 

953 

152.8 

168.8 

160.1 

87.0 

20 

80 

211 

847 

155.0 

165.6 

159.9 

86.9 

30 

70 

317 

741 

156.8 

163.0 

159.8 

86.8 

40 

60 

423 

635 

158.4 

161.2 

159.8 

86.7 

50 

50 

529 

529 

159.7 

159.7 

159.7 

86.7 

Table  1 1 .  Results  of  varying  extreme  radiance  thresholds  for  Solar  Village  north 
satellite  dish  sunlit  region  and  the  April  12,  2006  satellite  image. 


b.  Goddard  Space  Flight  Center  Building 

Figure  38  is  a  plot  of  the  sorted  spectral  radiance  of  usable  shadow  and 
sunlit  pixels  for  the  GSFC  building  using  the  January  16,  2009  image.  In  this  case,  the 
shadow  region  is  much  larger  than  the  sunlit  region  and  the  subsequent  slopes  are 
reversed  with  the  smaller  sunlit  region  (179  pixels)  being  steeper  than  the  shadow  region 
(751  pixels).  Otherwise,  the  shape  of  the  shadow  region  curve  is  similar  to  the  shadow 
region  curve  from  Solar  Village  with  the  highest  radiance  pixels  (32  of  748)  being 
significantly  brighter  than  the  rest  of  the  dataset.  Interestingly,  the  sunlit  region  curve  has 
a  different  shape  than  Solar  Village  with  the  brightest  pixels  being  much  brighter  than  the 
rest  of  the  sample  due  to  variations  in  surface  reflectance  from  two  distinct  surface  types; 
some  pavement  with  a  primarily  grassy  surface.  Otherwise,  the  impact  on  the  radiance 
difference  as  larger  percentages  of  the  extremes  from  the  total  sample  are  removed  is 
very  similar  to  the  Solar  Village  case  where  removing  10%  of  the  darkest  and  brightest 
pixels  has  the  largest  affect  on  the  average  radiances  (results  not  pictured). 

Although  this  suggests  that  10%  is  the  ideal  threshold  of  the  extremes,  the 
small  variation  as  more  of  the  extremes  are  rejected  leaves  open  the  possibility  that  larger 
portions  of  the  extremes  can  be  trimmed  without  significantly  changing  the  resulting 
AOD  values.  The  benefit  to  trimming  larger  percentages  of  the  extreme  radiances  should 
be  most  evident  when  calculating  AOD  over  heterogeneous  surfaces  or  when  overcoming 
issues  with  offset  correction  or  partial  shadow  pixel  rejection  thresholds.  The  results 
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from  the  GSFC  and  Solar  Village  tests  both  indicate  that  any  threshold  over  10%  should 
provide  representative  samples  of  the  shadow/sunlit  regions. 


+  Non-shadow 
-  Shadow 


Figure  38.  Spectral  radiance  of  usable  pixels  for  shadow  and  sunlit  regions  for  GSFC 

building  and  the  January  16,  2009  image. 


Figure  39  displays  the  results  from  running  FASTEC  on  the  GSFC  DSM 
for  the  October  26,  2006  image  with  different  thresholds  of  extreme  radiance  rejections. 
In  each  case,  both  extremes  were  trimmed  an  equal  percentage  for  both  the  shadow  and 
sunlit  regions.  This  graph  highlights  the  fact  that  computed  AOD  values  are  very 
consistent  for  each  of  the  nine  usable  shadows  in  this  image,  regardless  of  the  threshold 
setting.  In  this  scene,  3  of  the  9  AOD  values  are  consistently  negative  regardless  of  the 
threshold.  These  negative  values  occur  when  molecular  Rayleigh  optical  depth  is 
subtracted  from  total  optical  depth.  Otherwise,  computed  AOD  values  for  each  building 
vary  by  less  than  0.02  as  the  extreme  rejection  values  change.  These  facts  suggest  the 
threshold  has  little  impact  on  the  results,  however,  the  true  AOD  for  this  scene  is  0.059, 
well  below  the  recommended  threshold  of  0.10  (Vincent  2006).  In  addition,  the 
minimum  building  height  setting  happens  to  select  buildings  whose  shadows  fall  on 
largely  homogeneous  surfaces,  reducing  the  impact  of  these  variations  of  the  extreme 
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rejection  threshold.  Based  on  the  findings  from  this  section,  the  extreme  rejection 
thresholds  will  be  set  to  25%  for  the  whole  DSM  tile  tests  in  the  next  section.  This 
ensures  the  middle  50%  of  radiances  in  the  shadow/sunlit  regions  are  used  to  compute 
optical  depth  and  should  minimize  variations  in  AOD,  making  the  results  more  precise. 
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Figure  39.  Aerosol  optical  depth  using  different  extreme  radiance  rejection  percentages 
on  the  GSFC  DSM  file  and  the  October  26,  2006  image.  For  each  threshold,  the 
box  is  the  50th  percentile,  the  error  bars  are  the  25th  and  75th  percentile,  plus 
symbols  mark  the  maximum  and  minimum,  and  the  dotted  line  marks  ground 

truth  AOD. 

B.  WHOLE  SCENE  AOD  CALCULATION 
1.  Solar  Village 

The  first  full  scene  investigation  of  FASTEC  was  conducted  using  the  manually 
created  Solar  Village  elevation  file  and  14  satellite  images  (5  multispectral,  9 
panchromatic)  with  corresponding  ground  truth  AERONET  data.  This  scene  consisted  of 
simple  elevation  data  with  high  surface  reflectance  in  a  desert  region.  Figure  40  is  the 
manually  created  64  cm  horizontal  resolution  elevation  file  covering  an  approximately 
3600  square  meter  area  about  200  meters  northeast  of  the  AERONET  station  (located  at 
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N  24°  54’  25”  N,  E  46°  23’  49”).  The  surface  elevation  was  set  to  1  m,  the  rectangular 
building  elevation  was  set  to  5  m,  and  the  satellite  elevations  were  set  to  about  14  meters. 
Appendix  D  lists  the  configuration  settings  with  minimum  building  height  set  to  six  to 
ensure  AOD  calculations  were  only  performed  on  the  two  satellite  dishes,  not  the  small 
rectangular  building.  The  satellite  dishes  were  movable  and  their  orientation  was 
occasionally  different  than  pictured. 


Figure  40.  Manually  created  DSM  for  Solar  Village.  The  (1)  north  and  (2)  south  satellite 
dishes  were  estimated  to  be  13  meters  tall  and  the  building  was  estimated  to  be  4 

meters  tall. 

a.  Multispectral  Images 

Figure  41  displays  the  red,  green,  and  blue  (RGB)  channels  from  the  five 
multispectral  QuickBird  images  that  were  evaluated  using  FASTEC.  The  NIR  channel 
was  also  evaluated  (not  pictured).  Ground  truth  best  fit  curves  were  determined  with 
level  2.0  AERONET  daily  average  AOD  whenever  it  changed  by  less  than  0.02.  In  the 
case  of  the  April  30,  2006  image,  level  2.0  AERONET  data  stopped  at  0707  and  the 
image  was  collected  at  0755.  Level  1.0  data  showed  AOD  increasing  rapidly  during  the 
satellite  collection  time.  Surface  visibility  at  Riyadh  International  Airport,  located 
approximately  15  miles  southeast  of  the  test  location,  decreased  from  8  km  to  4  km 
between  0800  and  1000  (the  0900  observation  is  missing).  With  a  westerly  wind  of  8-10 
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knots  during  that  period,  it  would  take  between  one  and  two  hours  for  dust  to  advect  over 
Riyadh  from  the  region  where  the  AERONET  site  is  located.  Based  on  Koschmieder’s 
equation  (Horvath  1971),  these  AOD  values  are  consistent  with  observed  surface 
visibility  (see  Equations  (20)  and  (21)): 


Vis  = 


3.912 

°e 

5„ 


a.  = 


H 


(20) 

(21) 


where  5a  is  AOD,  ae  is  the  extinction  coefficient,  and  H  is  the  height  of  the  planetary 


boundary  layer.  For  this  day,  the  mean  of  the  two  surrounding  level  1.0  AERONET 
samples  was  used  instead  of  the  daily  average.  In  the  October  19,  2005  and  the 
November  30,  2007  images,  satellite  dish  orientations  were  noticeably  different  from  the 
elevation  file. 
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Figure  41.  Multi-spectral  imagery  for  Solar  Village  testing  on  (A)  October  19,  2005,  (B) 
April  12,  2006,  (C)  April  30,  2006,  (D)  May  26,  2007,  and  (E)  November  30, 
2007  with  2.4  meter  resolution  QuickBird  data. 


Basic  data  for  each  image  test  is  listed  in  Table  12.  Relatively  large 
shadows  (with  a  pixel  count  of  more  than  100)  were  identified  in  the  October  and 
November  images  while  the  April  and  May  images  had  relatively  small  shadows.  Offset 
correction  shifted  the  DSM  less  than  20  meters  in  any  direction  for  each  image.  In  the 
May  26  image  the  mean  sun  azimuth  and  elevation  information  in  the  metadata  file  that 
accompanied  the  image  was  incorrect  (150.2  and  44.9  respectively)  and  was  adjusted  by 
hand  based  on  historical  data  from  the  Naval  observatory. 


70 


Date 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Satellite 

Dish 

Shadow 

Pixel  Ct 

Non-shadow 

Pixel  Ct 

10/19/2005 

0752 

160.2 

53.4 

-17 

31 

1 

113 

643 

2 

101 

658 

4/12/2006 

0754 

135.4 

68.4 

3 

-2 

1 

27 

551 

2 

33 

588 

4/30/2006 

0755 

124.3 

73.5 

-8 

7 

1 

13 

476 

2 

11 

576 

5/26/2007 

0801 

105.2 

77.6 

-4 

-13 

1 

6 

382 

2 

9 

553 

11/30/2007 

0754 

164.9 

42.2 

-2 

29 

1 

107 

469 

2 

101 

560 

Table  12.  Date,  time,  sun  azimuth,  offset  correction,  and  pixel  counts  for  the  satellite 
dish  shadows  in  the  Solar  Village  multispectral  imagery. 


In  general,  AOD  values  using  FASTEC  were  quite  accurate  with  bias  less 
than  0.05  and  standard  deviation  less  than  0.08  (see  Figures  42-45  and  Table  13).  Using 
Equations  (20)  and  (21),  this  standard  deviation  equates  to  a  horizontal  visibility  variation 
of  less  than  1.5  km  when  AOD  is  0.5  or  greater  and  the  planetary  boundary  layer  height 
is  1  km.  The  largest  error  was  found  in  the  two  April  images.  On  April  30,  2006,  level 
1.0  data  showed  rapidly  increasing  AOD  to  over  1.0  around  the  time  the  image  was 
collected.  It  is  interesting  that  FASTEC  showed  an  unusually  large  decrease  in  AOD 
between  channels  when  comparing  April  30,  2006  with  the  other  days.  This  may  be  due 
to  the  size  of  the  aerosols  present  or  due  to  problems  with  the  shadow  method  when 
optical  depth  is  high.  Otherwise,  the  10  data  points  showed  a  slight  high  bias  for  the  blue 
channel  and  a  slight  low  bias  for  the  others  with  each  channel  error  standard  deviation 
less  than  0.08. 
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Figure  42.  Blue  channel  AOD  results  for  Solar  Village  satellite  dishes. 
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Figure  43.  Green  channel  AOD  results  for  Solar  Village  satellite  dishes. 
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Figure  44.  Red  channel  AOD  results  for  Solar  Village  satellite  dishes. 
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Figure  45.  Near-Infrared  channel  AOD  results  for  Solar  Village  satellite  dishes. 
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Channel 

Bias 

Standard 

Deviation 

Blue 

0.043 

0.078 

Green 

-0.006 

0.078 

Red 

-0.008 

0.079 

NIR 

-0.003 

0.078 

Table  13.  Multispectral  channel-specific  error  bias  and  error  standard  deviation 

Solar  Village  shadows. 


b.  Panchromatic  Images 

Nine  panchromatic  images  over  Solar  Village  were  evaluated  with 
FASTEC  (see  Figure  46).  The  satellite  dishes  are  pointing  in  different  directions  in 
several  of  the  images,  most  notably  (F)  and  (I).  This  should  affect  the  accuracy  of  the 
shadow  found  using  the  hand-generated  Solar  Village  DSM  with  the  dish  orientation 
pointing  almost  straight  up.  Ground  truth  best  fit  curves  were  based  on  the  daily  average 
AERONET  value  for  each  image  except  for  the  April  30,  2006  image.  As  discussed  in 
the  previous  section,  level  1.0  AERONET  data  showed  AOD  increasing  rapidly  from 
about  0.7  to  1.0  around  the  time  of  the  satellite  image.  For  every  other  image,  the  AOD 
changed  slightly  with  standard  deviations  less  than  0.02  for  the  whole  day. 
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Figure  46.  Panchromatic  imagery  for  Solar  Village  testing  on  (A)  October  19,  2005,  (B) 
April  12,  2006,  (C)  April  30,  2006,  (D)  May  26,  2007,  (E)  November  30,  2007, 
(F)  June  26,  2008,  (G)  October  11,  2008,  (H)  October  14,  2008,  (I)  December  14, 

2008. 

As  in  the  multispectral  images,  FASTEC  computed  AOD  for  both  satellite 
dishes  in  each  image.  Basic  data  on  each  of  the  scenes  is  provided  in  Table  14.  Four 
images,  collected  between  April  and  June,  had  relatively  small  shadows  with  fewer  than 
33  pixels  while  the  others  identified  between  80  and  200  pixels.  The  horizontal 
resolution  of  WorldView-1  data  (50  cm)  was  slightly  higher  than  DSM  resolution  (64 
cm)  while  the  QuickBird  imagery  had  the  same  resolution  as  the  DSM.  As  in  the 
multispectral  case,  the  May  26,  2007  sun  azimuth  and  elevation  was  incorrectly  identified 
in  the  metadata  file  and  was  manually  corrected  to  105.2°  and  77.6°  respectively. 
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Date 

Sat 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Satellite 

Dish 

Shadow 

Pixel  Ct 

Non-shadow 

Pixel  Ct 

10/19/2005 

QB 

0752 

160.2 

53.4 

-14 

31 

1 

113 

643 

2 

101 

658 

4/12/2006 

QB 

0754 

135.4 

68.4 

4 

-1 

1 

27 

551 

2 

33 

588 

4/30/2006 

QB 

0755 

124.6 

73.5 

-38 

14 

1 

13 

477 

2 

11 

577 

5/26/2007 

QB 

0801 

105.2 

77.6 

-13 

-2 

1 

6 

382 

2 

9 

553 

11/30/2007 

QB 

0754 

164.9 

42.2 

4 

30 

1 

108 

469 

2 

101 

560 

6/26/2008 

WV1 

0736 

90.3 

71.5 

2 

10 

1 

18 

446 

2 

20 

475 

10/11/2008 

WV1 

0743 

154.6 

55.1 

-7 

17 

1 

94 

628 

2 

82 

635 

10/14/2008 

QB 

0750 

158.6 

54.9 

-25 

-4 

1 

93 

614 

2 

83 

646 

12/14/2008 

WV1 

0746 

161.3 

39.7 

4 

22 

1 

199 

505 

2 

167 

488 

Table  14.  Date,  satellite,  time,  sun  azimuth  and  elevation,  offset  correction,  and 

pixel  counts  for  the  satellite  dish  shadows  in  the  Solar  Village  panchromatic 

imagery. 


Figure  46  is  a  plot  comparing  the  AOD  values  found  using  the  satellite 
images  with  AERONET.  The  results  showed  a  low  bias  of  0.136  and  an  error  standard 
deviation  of  0.071  for  all  18  AOD  calculations.  If  the  suspect  April  30,  2006  case  is 
removed  from  the  sample,  the  low  bias  becomes  0.115  and  the  error  standard  deviation 
lowers  to  0.042.  Figure  48  displays  the  bias  corrected  results.  Although  the  sample  is 
rather  small  and  the  bias  and  accuracy  are  found  using  the  one  scene,  this  figure  shows 
good  agreement  between  ground  truth  and  FASTEC. 
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Figure  47.  Panchromatic  channel  AOD  results  for  Solar  Village  satellite  dishes. 
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Figure  48.  Panchromatic  channel  bias-adjusted  AOD  results  for  Solar  Village  satellite 

dishes. 


77 


2. 


Goddard  Space  Flight  Center 


The  first  full  scene  investigation  of  FASTEC  using  a  BuckEye  DSM  was 
conducted  with  the  approximately  5  km  x  5  km  Goddard  Space  Flight  Center  (GSFC) 
file.  The  region  provides  a  mixture  of  urban  and  rural  terrain  with  large  buildings, 
asphalt,  and  forested  areas  (see  Figure  49).  Appendix  E  lists  the  configuration  values 
where  the  minimum  building  height  and  minimum  pixel  count  were  set  up  to  eliminate 
tree  shadows.  Using  these  settings  prevented  AOD  calculations  on  all  but  nine  buildings 
from  the  DSM,  including  the  Goddard  Space  Flight  Center  building  where  the 
AERONET  station  was  located  and  where  single  shadow  testing  in  Chapter  IV,  Section  B 
was  performed  (see  upper  right  hand  corner  of  Figure  49). 


Figure  49.  BuckEye  DSM  tile  over  the  Goddard  Space  Flight  Center  with  zoomed  in 
view  of  three  buildings  identified  as  shadow  generators  by  FASTEC.  Green  dots 
are  shadow  regions  with  sun  to  SE.  Red  dot  marks  location  of  AERONET  site. 
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a.  Multispectral  Images 

Two  multispectral  images  that  covered  the  GSFC  region  and  coincided 
with  AERONET  data  were  evaluated  with  FASTEC.  Figure  50  shows  the  RGB  images 
for  the  zoomed  in  region  of  Figure  49.  The  fourth  channel,  at  the  NIR  wavelength,  was 
also  evaluated  (not  pictured).  Ground  truth  best  fit  curves  were  based  on  the  daily 
average  AERONET  value  for  each  image  because  AOD  changed  less  than  0.02.  Portions 
of  the  AERONET  data  was  missing  from  the  cloud  screened  data  on  February  8,  2009.  It 
is  easy  to  identify  the  building  shadows  in  the  October  26,  2006  image  on  the  left  while 
the  shadows  in  the  February  8,  2009  image  on  the  right  are  covered  with  thick  cumulus 
clouds. 


Figure  50.  Examples  of  multispectral  imagery  for  GSFC  BuckEye  file  testing  from 
October  26,  2006  (left)  and  February  8,  2009  (right).  Images  are  approximately 
equivalent  to  the  zoomed  in  region  in  Figure  49. 

On  February  8,  2009,  every  usable  target  shadow  is  obscured  by  cloud. 
For  4  of  the  6  target  shadows,  FASTEC  found  shadow  radiance  exceeded  sunlit  radiance 
and  AOD  was  not  calculated.  The  other  target  shadows  calculated  unrealistic  AOD 
values,  between  0.43  and  1.13  with  ground  truth  of  0.07  for  the  red  channel  for  instance. 
For  this  reason,  the  February  8,  2009  image  was  not  included  in  the  results. 

On  October  26,  2006  at  1621Z,  the  sun  azimuth  and  elevation  were  170.7° 
and  38.2°  respectively.  Offset  correction  shifted  the  DSM  file  eight  meters  east  to  align  it 
with  the  satellite  pixels.  For  the  eight  usable  shadows,  the  average  shadow  region  was 
273  pixels  while  the  average  sunlit  region  was  296  pixels.  The  50th  percentile  value  of 
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AOD  for  the  blue  channel  was  slightly  above  ground  truth  while  the  50th  percentile  value 
was  below  ground  truth  for  the  other  channels  (see  Figure  51).  For  the  red  channel,  the 
calculated  AOD  was  less  than  zero  for  each  of  the  eight  shadows  and  the  50th  percentile 
was  0.091  below  ground  truth  (see  Table  15).  For  each  channel,  the  standard  deviation 
was  0.078  or  better.  For  the  blue  channel,  the  bias  was  noticeably  above  ground  truth  due 
to  one  shadow,  which  found  AOD  to  be  0.295.  Removing  that  shadow  from  the  dataset 
brought  the  bias  to  within  0.003  of  ground  truth  and  reduced  the  standard  deviation  to 
0.051. 


♦  Blue 

■  Green 
▲  Red 

■  NIR 


Figure  5 1 .  50th  percentile  value  of  AOD  vs.  AERONET  value  for  GSFC  multispectral 
channels  from  October  26,  2006  QuickBird  image.  Error  bars  on  the  y-axis 
reflect  the  reported  accuracy  of  AERONET  data  (0.02)  and  error  bars  on  x-axis 
are  one  standard  deviation  of  the  values  reported  on  that  day  and  scene.  n=8. 
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Channel 

Bias 

Standard 

Deviation 

Blue 

0.027 

0.078 

Green 

-0.055 

0.023 

Red 

-0.091 

0.039 

NIR 

-0.018 

0.020 

Table  15.  Multispectral  channel-specific  error  bias  and  standard  deviation  for 

October  26,  2006  GSFC  shadows. 

Interesting  results  appear  when  plotting  all  other  channels  against  the  blue 
channel  (see  Figure  52).  It  is  expected  that  AOD  would  decrease  with  increasing 
wavelength,  however,  FASTEC  AOD  values  decreased  from  blue  to  green  to  NIR  to  red. 
The  two  exceptions  were  the  cases  where  blue  channel  AOD  was  lowest.  In  the  case 
where  the  green  and  blue  channel  had  the  lowest  AOD  values,  the  sunlit  region  included 
water.  In  the  other  case,  green  was  lower  than  NIR  but  higher  than  red  and  the  sunlit 
region  was  composed  of  pavement,  grass,  and  cars. 

The  maximum  and  minimum  blue  channel  AOD  cases  are  also  noticeably 
different  than  the  other  six  cases.  These  extremes  were  primarily  due  to  different  surface 
types  dominating  the  sunlit  region.  In  the  maximum  blue  channel  AOD  case,  the  sunlit 
region  surface  consisted  primarily  of  grass  with  two  sidewalks  and  a  dirt  area  present. 
This  resulted  in  sunlit  radiance  being  significantly  lower  than  the  six  pavement  cases 
(36.5  RU  vs.  44-50  RU).  In  the  minimum  blue  channel  AOD  case,  the  sunlit  region 
surface  was  primarily  composed  of  sidewalk  and  water  with  some  grass  as  well.  This 
resulted  in  sunlit  radiance  being  significantly  higher  than  the  six  pavement  cases  (58  RU 
vs.  44-50  RU).  In  both  the  blue  channel  extremes,  the  green  and  NIR  AOD  was 
consistent  with  the  other  shadows.  For  the  maximum  blue  channel  AOD,  the  red  channel 
AOD  was  much  lower  than  the  other  shadows  because  its  sunlit  radiance  was 
significantly  lower  (19  RU  vs.  25-30  RU).  Overall,  shadow  radiance  varied  by  about  4 
RU  for  each  channel  except  NIR  where  it  varied  by  8  RU. 


81 


0.1  n 


0.05 


o 

o 

< 

0 

c 

c 

03 

-C 

o 

1— 

0 


-0.05 


-0.1 


♦ 

A 


A 


♦ 

▲ 


♦  Blue-Green 
■  Blue-red 
a  Blue-NIR 


-0.15 

0 


0.05  0.1  0.15  0.2  0.25  0.3  0.35 

Blue  channel  AOD 


Figure  52.  Plot  of  blue  channel  AOD  vs.  the  other  channels  for  the  October  26,  2006 

GSFC  image. 


b.  Panchromatic  Images 

Figure  53  shows  the  four  panchromatic  images  that  were  evaluated  with 
FASTEC.  Surface  obscuring  clouds  were  not  apparent  for  any  of  the  valid  target 
shadows  in  the  four  satellite  images  although  thin  cirrus  may  have  been  present.  Ground 
truth  was  determined  using  best  fit  curves  based  on  the  daily  average  AERONET  value 
for  each  image  because  AOD  changed  by  less  than  0.02  for  each  wavelength.  As 
explained  in  the  multispectral  section,  cloud  coverage  over  every  target  shadow  in  the 
February  8,  2009  led  to  unusable  results  and  the  image  was  discarded. 
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Figure  53.  Examples  of  panchromatic  imagery  for  GSFC  BuckEye  file  testing  from  (A) 
October  26,  2006,  (B)  March  5,  2008,  (C)  May  25,  2008,  and  (D)  January  16, 
2009.  Images  are  approximately  equivalent  to  the  zoomed  in  region  in  Figure  49. 


In  each  image,  FASTEC  computed  AOD  for  anywhere  from  three  to  eight 
shadows,  depending  on  the  overlap  between  the  DSM  file  and  the  satellite  image.  Basic 
data  is  listed  in  Table  16  where  the  top  image  was  from  QuickBird,  the  bottom  three  from 
WorldView-1.  Generally,  shadow  size  was  largest  for  the  January  16,  2009  image  where 
the  sun  elevation  was  lowest.  Offset  correction  shifted  the  DSM  between  about  3  and  14 
meters.  The  sun  azimuth  and  elevation  provided  in  the  metadata  files  appeared  accurate 


when  subjectively  compared  with  the  satellite  images. 


Date 

Time 

Sun 

Azimuth 

Sun 

Elevation 

Row 

Offset 

Column 

Offset 

Shadow 

Count 

Avg.  Pix 

Shadow 

el  Count 

Non-shadow 

10/26/2006 

1621 

170.7 

38.2 

2 

9 

8 

273 

296 

3/5/2008 

1543 

148.1 

40.4 

-2 

3 

7 

215 

308 

5/25/2008 

1606 

140.8 

68.4 

-1 

-13 

4 

79 

332 

1/16/2009 

1600 

159.6 

27.7 

-7 

5 

3 

560 

365 

Table  16.  Date,  time,  sun  azimuth  and  elevation,  offset  correction,  shadow  count, 

and  average  pixel  counts  for  the  satellite  dish  shadows  in  the  Solar  Village 

panchromatic  imagery. 
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Figure  54  displays  the  results  for  the  four  test  cases.  In  each  case,  AOD 
values  were  less  than  0.1,  with  FASTEC  finding  the  lowest  values  for  the  May  25,  2008 
image.  The  March  5  and  January  16  images  showed  the  50th  percentile  AOD  value 
slightly  above  ground  truth  while  the  other  two  images  were  below  ground  truth.  The 
standard  deviation  of  the  AOD  values  of  the  May  25  image  was  slightly  less  than  0.03 
while  the  other  images  were  less  than  0.02.  Overall,  the  four  test  cases  suggested  a  low 
bias  of  approximately  0.04.  In  the  May  25,  2008  case,  the  negative  AOD  values  occurred 
after  removing  the  Rayleigh  contribution. 


■  20061026 

♦  20080305 
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♦  20090116 


Figure  54.  50th  percentile  value  of  AOD  vs.  AERONET  value  for  GSFC  panchromatic 
images.  Error  bars  on  the  y-axis  reflect  the  reported  accuracy  of  AERONET  data 
(0.02)  and  error  bars  on  x-axis  are  one  standard  deviation  of  the  values  reported 

for  that  day  and  scene. 


3.  Other  Location  Tests 


a.  Tucson 

Tucson  was  chosen  to  test  FASTEC  using  BuckEye  DSM  data  in  a  desert 


region.  Figure  55  shows  the  DSM  file  used  for  this  location  with  the  zoomed  in  region 
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showing  the  area  surrounding  the  AERONET  station.  The  configuration  values  were  the 
same  as  those  used  for  GSFC  testing  (see  Appendix  E).  The  large  bright  region  in  the 
southwest  corner  of  the  DSM  is  caused  by  mountainous  terrain. 


Figure  55.  DSM  file  for  Tucson  testing.  The  inset  is  the  area  around  the  AERONET 
station,  the  bright  area  in  the  southwest  corner  is  caused  by  mountains. 

Only  one  available  Tucson  image  corresponded  with  AERONET  data. 
The  image  for  that  day,  November  1,  2005,  was  comprised  of  the  panchromatic  and  blue 
channels  from  the  IKONOS  satellite  (see  Figure  56).  The  four  meter  resolution  blue 
channel  is  much  fuzzier  than  the  one  meter  resolution  panchromatic  channel.  Ground 
truth  was  determined  using  best  fit  curves  of  the  daily  average  AERONET  AOD  because 
it  varied  by  less  than  0.02  for  each  wavelength. 
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Figure  56.  November  1,  2005  1815Z  imagery  for  Tucson.  Image  (A)  is  the  four-meter 
resolution  blue  channel  and  (B)  is  the  one-meter  resolution  panchromatic  channel. 

The  area  in  each  image  roughly  matches  the  zoomed  in  area  in  Figure  55. 

FASTEC  computed  AOD  for  both  channels.  The  sun  azimuth  and 
elevation  were  163°  and  41.6°  respectively.  Offset  correction  adjusted  the  DSM  13  rows 
and  24  columns  for  the  blue  channel  and  14  rows  and  24  columns  for  the  panchromatic 
channel.  Calculations  were  performed  on  14  shadows  in  the  blue  channel  with  an 
average  shadow  region  size  of  257.3  pixels  and  average  sunlit  region  size  of  256.4  pixels 
while  3  shadows  in  the  panchromatic  channel  had  an  average  size  of  199.3  shadow  pixels 
and  214.7  sunlit  pixels.  The  disparity  in  the  number  of  targets  between  the  two  channels 
was  because  the  areal  extent  of  the  panchromatic  channel  was  much  smaller  than  the  blue 
channel. 

The  blue  channel  AOD  values  had  a  significant  high  bias  of  over  0.28  (see 
Figure  57).  The  standard  deviation  of  AOD  values  for  the  blue  channel  was  0.055, 
suggesting  the  method  found  values  in  good  agreement  with  each  other  despite  being 
very  high.  It  appears  the  high  bias  was  caused  by  the  poor  resolution  collecting  data  over 
a  larger  area  for  each  pixel  in  the  satellite  image.  This  combined  with  the  weighting 
method  used  to  assign  radiance  values  to  the  DSM  pixels  smoothed  the  pixel  radiances, 
decreasing  the  radiance  difference  between  shadow/sunlit  regions,  and  causing  AOD  to 
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be  too  high.  The  complex  city  scene  may  amplify  this  problem  and  IKONOS 
multispectral  imagery  may  work  properly  in  scenes  with  more  homogeneous  surfaces. 
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Blue  channel  error  bin  (computed  -  truth) 

Figure  57.  Histogram  of  blue  channel  AOD  computation  error  for  November  1,  2005 

image.  Ground  truth  AOD  was  0.046. 

The  panchromatic  channel  expected  AOD  had  a  slight  low  bias  of  0.029 
relative  to  ground  truth  of  0.026,  with  each  computed  AOD  value  being  slightly  less  than 
zero  (-0.0045,  -0.0021,  -0.0008).  The  standard  deviation  of  AOD  values  for  the 
panchromatic  channel  was  0.001,  showing  strong  agreement.  The  findings  for  this 
channel  support  the  belief  that  the  comparatively  poor  resolution  of  the  IKONOS 
multispectral  sensor  may  be  insufficient  for  this  scene. 

b.  Baltimore  and  Houston 

Baltimore  and  Houston  were  included  to  test  FASTEC  in  very  complex 
settings  with  numerous  shadow  generators  located  in  close  proximity  to  each  other  (see 
Figure  58  and  59  respectively).  Of  the  three  Houston  images  tested,  none  had  coinciding 
ground  truth  AERONET  data  while  14  of  the  15  Baltimore  images  coincided  with 
AERONET  data.  Due  to  problems  with  FASTEC  in  these  complex  settings,  comparisons 
with  ground  truth  were  not  perfonned. 
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Figure  58.  Baltimore  DSM.  The  zoomed  in  region  shows  good  shadow  generators  near 

the  Maryland  Science  Center  AERONET  station. 


Figure  59.  Houston  DSM.  The  zoomed  in  region  shows  good  shadow  generators  near 

the  University  of  Houston  AERONET  station. 
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The  Houston  DSM  included  the  downtown  area  with  a  complicated 
shadow  patterns  (see  Figure  60).  Although  FASTEC  appeared  to  correctly  identify 
shadow  regions,  the  offset  correction  algorithm  failed  to  correctly  align  the  DSM  and 
satellite  file.  This  led  to  widely  varying,  inaccurate  AOD  values.  For  example,  the 
January  25,  2004  image  found  24  targets  for  AOD  calculation.  3  of  the  24  had  shadow 
radiance  greater  than  sunlit  radiance  and  AOD  was  not  calculated.  The  other  2 1  targets 
had  an  average  AOD  in  the  blue  channel  of  0.338  and  a  standard  deviation  of  0.449. 
Although  ground  truth  was  not  available,  the  comparatively  large  standard  deviation 
suggested  the  results  were  invalid.  When  the  locations  of  the  shadows  in  the  output  file 
were  compared  with  the  satellite  image,  it  was  obvious  offset  correction  failed  to 
correctly  align  the  satellite  and  DSM.  Attempts  to  align  them  by  hand  were  not 
successful  either.  The  other  Houston  images  suffered  from  poor  offset  correction  as  well. 


Figure  60.  Multispectral  QuickBird  image  of  Houston  downtown  area  from  January  25, 
2004  at  1647Z.  Note  the  complexity  of  shadows  cast  by  the  tall  buildings. 
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The  Baltimore  DSM  included  the  downtown  area  and  Chesapeake  Bay 
(see  Figure  61).  The  large  body  of  water  further  complicated  the  scene  and  contributed  to 
the  failure  of  offset  correction  and  AOD  values  to  be  wrong.  It  is  theorized  that  the  low 
radiance  over  water  offered  a  large  area  where  the  difference  between  shadow  and  sunlit 
radiance  was  quite  small,  contributing  to  the  failure  of  the  offset  correction  algorithm  to 
properly  adjust  the  DSM.  For  the  January  1,  2004  panchromatic  image,  the  average 
AOD  was  0.305  and  the  standard  deviation  was  0.345  suggesting  the  results  were  invalid. 
When  the  location  of  shadows  in  the  output  file  were  compared  with  the  actual  shadows 
in  the  satellite  image,  it  was  apparent  that  offset  correction  again  failed.  In  fact,  offset 
correction  failed  to  correctly  align  the  DSM  and  satellite  file  in  every  Baltimore  and 
Houston  image  although  it  worked  properly  in  images  from  Solar  Village,  GSFC,  and 
Tucson.  Several  attempts  to  manually  correct  for  the  offset  failed  to  adequately  align  the 
images.  As  a  result,  the  results  from  complex  urban  scenes  were  discarded  for  this  study. 


Figure  61.  Panchromatic  QuickBird  image  of  Baltimore  from  January  1,  2004  at  1540Z. 
Note  the  complexity  of  shadows  cast  by  the  tall  buildings  and  the  optically  dark 

Chesapeake  Bay. 
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4.  Summary  of  Whole  Scene  AOD  Calculations 


FASTEC  showed  promise  when  used  in  the  Solar  Village  and  GSFC  scenes.  In 
both  scenes,  small  bias  pointed  to  problems  with  accuracy  but  FASTEC  calculated  AOD 
values  that  were  consistent  with  each  other,  as  evidenced  by  standard  deviations  less  than 
0.08.  Results  from  the  GSFC  scenes  suggest  that  bias  correction  may  work  best  when  the 
surface  type  is  known.  The  test  of  the  IKONOS  image  in  Tucson  suggested  the  lower 
resolution  multispectral  sensor  may  be  problematic,  at  least  in  urban  settings,  but  the 
panchromatic  channel  appeared  to  work  properly.  The  complex  scenes  in  Houston  and 
Baltimore  highlighted  problems  with  the  offset  correction  procedure  and  identified  the 
need  to  improve  upon  that  procedure  or  use  higher  quality  orthorectified  satellite  imagery 
before  testing  FASTEC  in  such  complex  urban  settings. 

C.  DSM  RESOLUTION  ADJUSTMENT 

The  final  test  of  FASTEC  was  designed  to  determine  the  minimum  usable 
resolution  of  elevation  data  by  reducing  the  horizontal  resolution  of  the  DSM.  The 
horizontal  resolution  of  the  DSM  was  reduced  using  two  methods.  Method  one  sampled 
a  subset  of  the  full  resolution  DSM  data  while  method  two  averaged  the  full  resolution 
DSM  data  (see  Figure  62).  The  first  method  is  designed  to  simulate  what  the  DSM  might 
look  like  if  a  LIDAR  sampled  the  surface  with  half  the  frequency  of  BuckEye  while  the 
second  method  was  designed  to  simulate  what  the  results  might  look  like  using  stereo 
pairs  DSM  data. 


91 


A 


15 

30 

30 

30 

15 

30 

30 

30 

15 

15 

15 

30 

15 

15 

15 

15 

15 

30 

15 

15 

22.5 

30 

15 

18.75 

Figure  62.  Example  of  resolution  adjustment  methods.  Box  A  represents  a  full  resolution 
DSM,  box  B  and  C  halve  the  horizontal  resolution  using  method  one  (sampling) 
and  two  (smoothing)  respectively.  Gray  shading  is  by  elevation. 


1.  Resolution  Reduction  for  Solar  Village  Scene 

Figure  63  shows  examples  of  reduced  resolution  DSM  images  using  the  sampling 
and  smoothing  method  with  resolution  reduced  to  1.28  and  1.92  meters.  The  most 
notable  difference  between  the  two  methods  is  the  edges  of  the  objects  are  less  sharp 
when  the  smoothing  method  is  applied.  The  most  notable  difference  between  the 
resolution  changes  with  either  method  is  the  larger  area  of  the  full  DSM  that  is  visible  in 
the  30  pixel  by  30  pixel  region.  Otherwise,  resolution  reduction  using  sampling  makes 
the  edges  appear  sharper  while  smoothing  makes  the  edges  less  sharp.  The  configuration 
values  for  this  test  are  listed  in  Appendix  D.  Resolution  reductions  beyond  1.92  meters 
lowered  the  pixel  count  below  the  user-defined  minimum  of  5. 
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Figure  63.  Examples  of  top  left  30  x  30  pixel  region  in  Solar  Village  DSM  after 

resolution  adjustment  methods  applied. 


Table  17  displays  the  results  of  resolution  reduction  testing  using  the  sampling 
method  with  the  Solar  Village  elevation  file  and  the  multispectral  satellite  channels.  As 
elevation  file  resolution  is  reduced,  AOD  values  decrease  with  the  overall  decrease 
typically  0.05  for  the  blue  channel  and  about  0.03  for  the  green,  red,  and  NIR  channels. 
These  decreases  are  primarily  due  to  decreasing  shadow  radiance.  For  example,  blue 
channel  radiance  values  decrease  by  about  1  RU  for  each  unit  of  resolution  decrease  on 
dish  1  and  3  RU  for  each  resolution  change  on  dish  2.  These  cause  an  increase  the 
radiance  difference  and  decrease  in  computed  AOD  values.  The  fact  that  the  changes  are 
consistent  for  each  shadow  suggests  lower  resolution  elevation  data  may  still  be  usable, 
with  only  the  bias  needing  adjustment. 
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Date 

Dish 

Channel 

AOD  w/  Resolution 

0.64 

1.28 

1.92 

Truth 

10192005 

1 

Blue 

0.301 

0.296 

0.256 

0.372 

Green 

0.241 

0.232 

0.215 

0.347 

Red 

0.223 

0.214 

0.195 

0.317 

NIR 

0.220 

0.212 

0.198 

0.297 

2 

Blue 

0.293 

0.268 

X 

0.372 

Green 

0.253 

0.232 

X 

0.347 

Red 

0.254 

0.231 

X 

0.317 

NIR 

0.254 

0.231 

X 

0.297 

11302007 

1 

Blue 

0.263 

0.247 

0.213 

0.168 

Green 

0.202 

0.181 

0.164 

0.149 

Red 

0.171 

0.155 

0.154 

0.124 

NIR 

0.157 

0.140 

0.144 

0.104 

2 

Blue 

0.266 

0.236 

X 

0.168 

Green 

0.217 

0.193 

X 

0.149 

Red 

0.196 

0.165 

X 

0.124 

NIR 

0.181 

0.156 

X 

0.104 

Table  17.  Results  of  Solar  Village  resolution  reduction  using  sampling  approach  and 
the  multispectral  channels.  x=insufficient  pixels  for  calculation. 


Table  18  displays  the  results  of  resolution  reduction  testing  using  the  sampling 
method  with  the  Solar  Village  elevation  file  and  the  panchromatic  satellite  channel.  With 
resolution  of  1 .92  meters,  AOD  was  not  calculated  on  the  southern  satellite  dish  (Dish  2) 
and  on  the  north  dish  on  October  11,  2008.  Otherwise,  AOD  varied  less  than  0.01  for 
every  image  except  the  November  30,  2007  image  where  it  varied  about  0.023.  Variation 
this  small  suggests  FASTEC  works  with  elevation  files  of  up  to  1 .92  meter  resolution  and 
possibly  beyond  since  FASTEC  stopped  calculating  AOD  before  the  values  changed 
significantly. 

Results  from  the  smoothing  method  (not  pictured)  were  very  similar  to  the 
sampling  method  with  one  exception.  On  October  14,  2008  with  resolution  at  1.92 
meters,  the  Dish  1  AOD  value  was  0.777  because  shadow  radiance  increased  from 
around  55  RU  to  121.4  RU  with  only  3  pixels.  Otherwise,  as  with  the  multispectral  case, 
the  AOD  value  changes  appear  consistent  for  each  shadow,  suggesting  using  lower 
resolution  data  may  only  affect  the  bias  of  the  results. 
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Date 

Dish 

AOD  w/  Resolution 

AERONET 

AOD 

0.64 

1.28 

1.92 

10112008 

1 

0.095 

0.093 

X 

0.267 

2 

0.096 

0.091 

X 

10142008 

1 

0.113 

0.109 

0.108 

0.245 

2 

0.115 

0.107 

X 

10192005 

1 

0.175 

0.167 

0.165 

0.319 

2 

0.183 

0.174 

X 

11302007 

1 

0.088 

0.084 

0.066 

0.124 

2 

0.087 

0.087 

X 

12142008 

1 

0.091 

0.090 

0.088 

0.230 

2 

0.104 

0.097 

X 

Table  18.  Results  of  Solar  Village  resolution  reduction  testing  using  sampling 

approach  and  the  panchromatic  channel.  x=insufficient  pixels  for  calculation. 


2.  Resolution  Reduction  for  GSFC  Scene 

Figure  64  shows  examples  of  reduced  resolution  DSM  images  using  the  sampling 
and  smoothing  method  with  resolution  reduced  to  2  and  3  meters.  Although  more 
difficult  to  discern  than  in  the  Solar  Village  image  in  Figure  63,  the  edges  of  the  objects 
are  less  sharp  when  the  smoothing  method  is  applied.  This  is  most  noticeable  with  the 
road  running  north-south  and  the  trees  south  of  the  building.  The  configuration  file  was 
set  to  the  values  listed  in  Appendix  E  with  the  exception  of  the  minimum  pixel  count 
variable  that  was  reduced  to  from  100  to  6.  This  change  was  made  to  allow  the  resolution 
reduction  algorithm  to  work  with  the  expected  lower  pixel  counts.  In  cases  where  this 
change  caused  a  tree  shadow  to  have  AOD  computed,  the  tree  shadow  was  removed  from 
the  sample. 
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Figure  64.  Examples  of  top  left  201  x  201  pixel  region  in  GSFC  DSM  after  resolution 

adjustment  methods  applied. 

Table  19  displays  the  results  of  resolution  reduction  testing  on  the  October  26, 
2006  GSFC  multispectral  image.  The  50th  percentile  AOD  values  tend  to  decrease  with 
reduced  elevation  file  resolution,  although  the  change  is  on  the  order  of  0.02  for  all 
channels.  These  decreases  are  primarily  caused  by  decreasing  shadow  radiance  that 
increases  the  radiance  difference  between  the  regions.  As  in  the  Solar  Village  test,  the 
fact  that  the  changes  are  consistent  for  each  shadow  suggests  lower  resolution  elevation 
data  may  still  be  usable,  with  only  the  bias  needing  adjustment.  This  test  suggests  4 
meter  horizontal  resolution  elevation  files  may  be  adequate  for  computing  AOD, 
although  lower  resolution  may  work  with  very  large  objects  in  simple  scenes. 
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Blue  AOD 

Green  AOD 

Red  AOD 

Number 

Range 

Range 

Range 

Range 

DSM 

of  calcs 

50% 

50% 

(75%-25%) 

50% 

(75%-25%) 

50% 

Full 

8 

0.121 

0.049 

0.031 

0.017 

-0.017 

0.027 

0.020 

0.024 

Sample  2 

6 

0.107 

0.016 

0.025 

0.011 

-0.017 

0.010 

0.017 

0.016 

Sample  3 

5 

0.106 

0.045 

0.021 

0.007 

-0.022 

0.006 

0.012 

0.011 

Sample  4 

5 

0.098 

0.004 

0.017 

0.015 

-0.037 

0.013 

0.003 

0.008 

Sample  5 

2 

0.084 

0.019 

0.007 

0.002 

-0.031 

0.004 

0.000 

0.007 

Smooth  2 

6 

0.110 

0.021 

0.026 

0.010 

-0.020 

0.023 

0.018 

0.010 

Smooth  3 

5 

0.095 

0.027 

0.016 

0.003 

-0.030 

0.004 

0.017 

0.010 

Smooth  4 

5 

0.088 

0.028 

0.014 

0.004 

-0.029 

0.009 

0.012 

0.003 

Smooth  5 

1 

0.095 

0.002 

-0.050 

-0.016 

Truth  0.104  0.083  0.058  0.040 


Table  19.  Results  of  October  26,  2006  GSFC  multispectral  resolution  reduction 

testing. 
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V.  CONCLUSIONS 


A.  SUMMARY 

FASTEC  worked  well  for  Solar  Village  and  the  Goddard  Space  Flight  Center 
testing  but  failed  to  work  properly  for  other  test  locations.  Shadow  selection  criteria 
testing  suggested  it  is  important  to  remove  partial  shadow  edge  pixels  from  the  sample. 
This  testing  also  suggested  it  is  somewhat  important  to  eliminate  pixels  that  are  close  to 
the  building  and  eliminate  a  portion  of  the  extremes  from  the  shadow  and  sunlit  regions. 
The  results  from  sky  dome  blockage  testing  were  inconclusive.  Whole  tile  testing  does 
suggest  it  is  useful  to  calculate  AOD  over  a  region  although  the  test  set-up  and  limited 
data  set  prevented  fully  assessing  the  accuracy  (bias  correction)  or  precision  (variation)  of 
FASTEC.  Finally,  the  DSM  resolution  test  showed  that  it  is  possible  to  use  lower 
resolution  elevation  files  and  still  achieve  good  results.  Further  testing  may  reveal  that 
the  elevation  file  resolution  must  only  be  capable  of  discerning  shadows  to  employ 
FASTEC. 

Strengths  of  FASTEC  include  the  speed  with  which  this  method  identifies  good 
targets,  the  objective  detennination  of  shadow/sunlit  regions,  the  easily  configurable 
nature  of  the  algorithm,  and  the  consistency  in  AOD  values  between  targets  in  the  same 
scene.  When  offset  correction  worked  properly,  the  technique  consistently  calculated 
AOD  values  within  0.05  for  targets  in  a  scene.  Elevation  data  resolution  appears  usable 
up  to  4  meters  depending  on  the  scene  being  analyzed.  This  technique  can  fill  the  high 
surface  reflectance  region  gap  in  remotely  sensed  AOD  databases.  Militarily,  FASTEC 
could  enable  in-scene  correction  in  intelligence  collection  and  possibly  influence  electro- 
optical  targeting  decisions. 

The  primary  weakness  of  this  study  is  the  failure  of  the  offset  correction 

algorithm  in  complex  urban  settings  such  as  Houston  and  Baltimore.  When  the  output  of 

FASTEC  has  high  standard  deviation,  the  user  should  locate  several  target  shadows  to 

determine  if  the  cause  of  the  high  standard  deviation  was  due  to  offset  correction  error, 

the  presence  of  clouds,  or  some  other  problem.  Additionally,  this  study  suggested  the 
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shadow  technique  has  bias,  which  may  vary  between  channels,  scenes,  or  surface  types. 
Although  only  a  few  images  were  used  in  this  study,  the  automated  approach  makes 
future,  larger  studies  using  different  sensors  feasible  for  shadow  analysis.  These  studies 
should  refine  the  bias  and  uncertainty  of  the  shadow  method. 

FASTEC  is  capable  of  supporting  operational  customers  or  scientists.  An 
intelligence  analyst  might  use  FASTEC  to  detennine  the  optical  depth  of  a  scene  and  use 
that  value  to  perform  in-scene  corrections  to  the  image.  A  climatologist  or  numerical 
modeler  might  establish  a  global  grid  of  shadow  generators  to  feed  a  global  aerosol 
database  while  the  ability  to  examine  detailed  shadow  characteristics  allows  researchers 
to  determine  different  shadow  characteristics. 

B.  FUTURE  RESEARCH 

Future  research  should  focus  on  two  areas:  testing  FASTEC  with  large  datasets  or 
different  sensors  and  improving  the  technique  itself.  When  it  comes  to  large  datasets,  test 
locations  should  include  AERONET  data  or  another  reliable,  accurate  means  of 
measuring  ground  truth.  Elevation  data  can  come  from  LIDAR  (as  was  used  in  this 
study),  stereo  pair  imagery,  or  other  sources  capable  of  providing  4  meter  horizontal 
resolution.  Using  large  datasets  may  enable  the  user  to  sort  results  by  surface  type, 
building  size,  or  sensor.  The  ability  to  define  small  sub-regions  within  the  DSM  should 
simplify  this  type  of  study.  If  the  user  analyzed  each  shadow  generator  in  one  full  scene 
run  and  created  configuration  files  around  the  desired  target  shadows,  the  output  could  be 
sorted  for  straightforward  statistical  analysis.  Similar  experiments  employing  different 
sensors  with  different  wavelength  bands  or  resolution  are  also  possible.  For  example, 
hyperspectral  images  provide  hundreds  of  usable  wavelength  bands  and  may  even  allow 
aerosol  size  characterization.  Additionally,  it  may  be  possible  to  use  lower  resolution 
satellites  exploiting  mountain  or  cloud  (if  cloud  top  heights  are  known)  shadows. 
Finally,  testing  FASTEC  in  complex  urban  settings  with  fully  orthorectified  satellite 
imagery  would  determine  the  feasibility  of  using  FASTEC  in  such  scenes. 

There  are  several  possible  ways  of  improving  FASTEC  to  include  eliminating 
dark  surfaces,  applying  a  cloud  mask,  identifying  isolated  buildings,  and  improving  offset 
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correction.  Optically  dark  surfaces  such  as  water  or  trees  pose  a  problem  for  FASTEC 
because  they  could  cause  offset  correction  to  fail  to  correctly  align  the  elevation  and 
satellite  files.  If  pixels  in  the  elevation  file  were  identified  by  surface  type,  these  regions 
could  be  avoided.  Similarly,  it  is  necessary  to  employ  the  shadow  technique  in  cloud-free 
regions.  Adding  a  cloud  mask  would  eliminate  cloud-obscured  targets  from  the 
calculation  and  improve  the  technique.  This  mask  should  be  applied  after  identifying 
shadow  pixels  but  before  performing  offset  correction.  A  robust  cloud  mask  should 
account  for  the  shadow  cast  by  the  cloud  and  block  off  a  standard  size  region  around  the 
cloud  to  eliminate  the  indirect  reflective  affects  of  the  cloud.  Also,  Vincent  (2006) 
suggested  selecting  isolated  shadows  for  AOD  calculation.  It  may  improve  the  technique 
to  eliminate  target  shadows  when  other  shadow  generators  within  a  user-defined  area  are 
present. 

Offset  correction  should  improve  when  working  with  fully  orthorectified  imagery. 
In  addition,  it  may  improve  with  the  addition  of  a  dark  surface  and  cloud  masks  but 
additional  changes  may  be  required.  First,  the  current  method  works  by  trying  to  align 
the  shadow  regions  in  the  DSM  with  low  radiance  pixels  in  the  satellite  file. 
Incorporating  a  dark  surface  mask  may  eliminate  this  issue  by  marking  these  low 
radiance  pixels  invalid  while  a  cloud  mask  would  prevent  high  radiance  cloud  pixels 
from  skewing  the  correction.  If  these  masks  fail  to  correct  the  problem,  an  alternative 
approach  is  to  adjust  the  DSM  position  to  maximize  the  radiance  difference  between 
shadow  and  sunlit  pixels.  This  approach  would  take  both  shadow  and  sunlit  radiance  into 
account  unlike  the  current  approach  that  just  examines  shadow  radiance.  Additionally,  it 
may  be  better  to  perform  offset  correction  on  each  individual  shadow  instead  of  the  entire 
elevation  file.  This  approach  would  increase  computation  time  but  make  each  AOD 
calculation  independent,  improving  the  likelihood  that  some,  or  most,  of  the  shadows  are 
properly  aligned.  Regardless,  correcting  for  the  offset  appears  to  be  the  biggest  obstacle 
to  implementing  this  algorithm  operationally. 
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APPENDIX  A:  FORTRAN  CONFIGURATION  VARIABLES 


A.  GENERAL  USE  VARIABLES  (CITY  VARS) 


Name 

Description 

Sat_path 

Path  to  satellite  file 

DSM_path 

Path  to  DSM  file 

op_path 

Path  to  output  directory 

Sat  name 

Base  input  satellite  filename 

DSM  name 

Base  input  DSM  filename.  Do  not  include  a  “.hdr”  file  descriptor 

Op  name 

Base  output  filename.  Makes  “AOD.txt”  and  “.geninfo.txt”  files 

DSM_sr/er, 

DSM_sc/ec 

Limits  shadow  search  to  sub  region  of  DSM  file.  If  set  outside 

the  bounds  of  the  DSM  file,  the  entire  DSM  file  is  explored. 

B.  SHADOW  REGION  VARIABLES  (SHADOW_VARS) 


Name 

Description 

DSM  min  bldg  hgt 

Minimum  height  of  shadow  generator 

DSMbldgdist 

Percentage  of  shadow  near  building  that  is  rejected 

DSM  min  bldg  dist 

Minimum  distance  of  valid  shadow  pixel  from  shadow  generator 

rlus/rlls 

Upper/lower  bounds  of  sorted  shadow  pixel  radiances 

max  ediff 

Maximum  range  of  elevation  differences  among  shadow  pixels 

C.  SUNLIT  REGION  VARIABLES  (SUNLIT  VARS) 


Name 

Description 

DSM  u  rgn 

Size  of  region  around  shadow  to  search  for  sunlit  pixels 
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nonshad  int 

Azimuth  interval  to  use  when  searching  for  sunlit  pixels 

nonshad  dint 

Distance  interval  on  each  azimuth  for  finding  sunlit  pixels 

nonshad  ediff 

Max  elevation  difference  of  sunlit  pixel  and  shadow 

rlun/rlln 

Upper/lower  bounds  of  sorted  sunlit  pixel  radiances 

D.  SHADOW/SUNLIT  REGION  VARIABLES  (SHADNONVARS) 


Name 

Description 

DSM  how  int 

Depth  of  shadow  edge  to  mark  invalid  due  to  edge  contamination 

minpix 

Minimum  count  of  shadow/sunlit  pixels  to  calculate  AOD 

E.  MISCELLANEOUS  VARIABLES  (MISC_VARS) 


Name 

Description 

whichAz 

Which  azimuth  (min/mean/max)  in  metadata  file  to  use 

whichEl 

Which  elevation  (min/mean/max)  in  metadata  file  to  use 

edgeadj 

Defunct 

off  cor 

Set  to  .true,  to  perform  automated  offset  correction 

roff,  c  off 

Used  to  manually  offset  DSM  row  and  column 

offshad 

Defunct 

offset 

Size  of  adjustment  region  for  the  DSM  file 

offchan 

Offset  correction  sat  channel  (set  to  1  if  value  exceeds  #  of  bands) 

chk 

Defunct 

RR 

Weighting  factor  for  nearest  neighbor  determination 

DSM  invalid 

Value  to  use  for  program  identified  invalid  pixels 

bigdatadump 

Set  to  .true,  to  dump  DSM  arrays  to  text  files  for  detailed  analysis 
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defunct 


parallel 


F.  VARIABLES  FOR  DOME  BLOCKAGE  CALCULATION  (DOME  VARS) 


Name 

Description 

error  check 

Set  to  .true,  to  perform  dome  blockage  calculation  on  each  pixel 

domesize 

Size  of  region  around  pixel  to  calculate  dome  blockage 

aziinc 

Azimuth  increment  for  dome  calculation 

distinc 

Distance  increment  for  dome  calculation 

maxdomes 

Maximum  allowed  dome  blockage  for  shadow  pixels 

maxdomen 

Maximum  allowed  dome  blockage  for  sunlit  pixels 

G.  AOD  CALCULATION  VARIABLES 


Name 

Description 

interval 

Scattering  angle  interval  in  MAR  calculation 

ssa 

Single-scatter  albedo  for  AOD  calculation 

g 

Henyey-Greenstein  constant 
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APPENDIX  B:  LIST  OF  BASELINE  USER-DEFINED  VARIABLES 
FOR  SHADOW  CHARACTERIZATION  TESTS  (SOLAR  VILLAGE) 


0.0 

DSM  min  bldg  hgt 

F 

off  cor 

1.0 

DSM_bldg_dist 

0 

r  off 

0.0 

DSM  min  bldg  dist 

0 

coff 

1.0 

Rlus 

100 

offshad 

0 

Rlls 

10 

offset 

1.5 

max  ediff 

1 

offchan 

15 

DSM  u  rgn 

T 

bigdatadump 

3 

nonshad  int 

T 

parallel 

0.99 

nonshad  dint 

T 

error  check 

1.0 

nonshad  ediff 

20 

domesize 

1.0 

rlun 

10 

aziinc 

0.0 

rlln 

2 

distinc 

1 

DSM  how  int 

1.0 

maxdomes 

15 

minpix 

1.0 

maxdomen 

2 

whichAz 

10 

interval 

3 

whichEl 

0.88 

ssa 

200 

edgeadj 

0.65 

g 
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APPENDIX  C:  LIST  OF  BASELINE  USER-DEFINED  VARIABLES 
FOR  SHADOW  CHARACTERIZATION  TESTS  (GSFC  BUILDING) 


15.0 

DSM  min  bldg  hgt 

F 

off  cor 

1.0 

DSM_bldg_dist 

0 

r  off 

0.0 

DSM  min  bldg  dist 

0 

coff 

1.0 

Rlus 

100 

offshad 

0 

Rlls 

10 

offset 

1.5 

max  ediff 

1 

offchan 

15 

DSM  u  rgn 

T 

bigdatadump 

3 

nonshad  int 

T 

parallel 

0.99 

nonshad  dint 

T 

error  check 

1.0 

nonshad  ediff 

20 

domesize 

1.0 

rlun 

10 

aziinc 

0.0 

rlln 

2 

distinc 

2 

DSM  how  int 

1.0 

maxdomes 

15 

minpix 

1.0 

maxdomen 

2 

whichAz 

10 

interval 

3 

whichEl 

0.88 

ssa 

200 

edgeadj 

0.65 

g 
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APPENDIX  D:  LIST  OF  USER-DEFINED  VARIABLES  FOR 
SOLAR  VILLAGE  AOD  CALCULATIONS 


6.0 

DSM  min  bldg  hgt 

T 

off  cor 

0.8 

DSM_bldg_dist 

0 

r  off 

2.0 

DSM  min  bldg  dist 

0 

coff 

0.75 

Rlus 

100 

offshad 

0.25 

Rlls 

40 

offset 

1.5 

max  ediff 

4 

offchan 

15 

DSM  u  rgn 

F 

bigdatadump 

3 

nonshad  int 

T 

parallel 

0.99 

nonshad  dint 

F 

error  check 

1.0 

nonshad  ediff 

20 

domesize 

0.75 

rlun 

10 

aziinc 

0.25 

rlln 

2 

distinc 

1 

DSM  how  int 

1.0 

maxdomes 

5 

minpix 

1.0 

maxdomen 

2 

whichAz 

10 

interval 

3 

whichEl 

0.88 

ssa 

200 

edgeadj 

0.65 

g 
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APPENDIX  E:  LIST  OF  USER-DEFINED  VARIABLES  FOR 
BUCKEYE  FILE  AOD  CALCULATIONS 


35.0 

DSM  min  bldg  hgt 

T 

off  cor 

0.5 

DSM_bldg_dist 

0 

r  off 

2.0 

DSM  min  bldg  dist 

0 

coff 

0.75 

Rlus 

100 

offshad 

0.25 

Rlls 

150 

offset 

1.5 

max  ediff 

4 

offchan 

15 

DSM  u  rgn 

T 

bigdatadump 

3 

nonshad  int 

T 

parallel 

0.99 

nonshad  dint 

F 

error  check 

1.0 

nonshad  ediff 

20 

domesize 

0.75 

rlun 

10 

aziinc 

0.25 

rlln 

2 

distinc 

2 

DSM  how  int 

1.0 

maxdomes 

100 

minpix 

1.0 

maxdomen 

2 

whichAz 

10 

interval 

3 

whichEl 

0.88 

ssa 

200 

edgeadj 

0.65 

g 
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